Blunt-nosed Leopard Lizard Data Analysis: Water Physiology Nuances

Savannah Weaver

Packages

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
if (!require("tidyverse")) install.packages("tidyverse") 
library("tidyverse") # for tidy/dplyr workflow of calculations
if (!require("lmerTest")) install.packages("lmerTest") 
library("lmerTest") # LMMs
if (!require("lme4")) install.packages("lme4") 
library("lme4") # LMMs
if (!require("Matrix")) install.packages("Matrix") 
library("Matrix") # LMMs
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # pretty colors
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("emmeans")) install.packages("emmeans")
library("emmeans") # marginal means & confints

#update.packages()

Set Colors & Themes

savs_ggplot_theme <- theme_classic() + 
                     theme(text = element_text(color = "black", 
                                               family = "sans", 
                                               size = 12),
                           axis.text = element_text(color = "black", 
                                                    family = "sans", 
                                                    size = 8),
                           legend.text = element_text(color = "black", 
                                                      family = "sans", 
                                                      size = 8),
                           legend.text.align = 0,
                           legend.position = "bottom"
                           )
  
# custom colors
Male <- brewer.pal(11, "RdYlBu")[c(10)]
Fem <- brewer.pal(11, "RdYlBu")[c(3)]
Fem_grav <- brewer.pal(11, "RdYlBu")[c(2)]
Fem_nongrav <- brewer.pal(11, "RdYlBu")[c(4)]
my_colors <- c(Fem, Male)
my_shapes <- c(19,15)
my_shapes_open <- c(21,22)

# corr colors
n3 <- brewer.pal(11, "RdYlBu")[c(9, 11, 3)]
n2 <- brewer.pal(9, "YlGnBu")[c(6, 8)]
n1 <- brewer.pal(9, "YlGnBu")[c(6)]


# labels for plots for repeat measures lizards only
my_labels <- c("April\nF\n(n=11)", "April\nM\n(n=22)",
               "May\nF\n(n=10)", "May\nM\n(n=17)",
               "July\nF\n(n=2)", "July\nM\n(n=12)")
# labels for plots for full dataset of all lizards
my_labels_full <- c("April\nF\n(n=26)", "April\nM\n(n=41)",
                     "May\nF\n(n=15)", "May\nM\n(n=22)",
                     "July\nF\n(n=3)", "July\nM\n(n=12)")
my_labels_full_F <- c("April\n(n=26)",
                     "May\n(n=15)",
                     "July\n(n=3)")
my_labels_full_M <- c("April\n(n=41)",
                     "May\n(n=22)",
                     "July\n(n=12)")

Background & Introduction

This data was collected on Blunt-nosed Leopard Lizards (Gambelia sila) during their active season in 2021. In this document, we will investigate the nuances and trends in the hydric physiology of these desert lizards.

We will do some additional data wrangling (see other ‘wrangling_’ files in this repository for more details on data preparation), compute some initial summary statistics, and investigate the following questions:

  • Are the different hydric physiology metrics we measured related to each other?
  • Did our attempt at supplemental hydration actually do anything?
  • How did hydric physiology change over time throughout the active season (April-July) of these lizards? This section composes the biggest chunk of this file.
  • Does hydric physiology relate to reproductive ability or effort in any way?

Load in Data

Main Dataframe

dat <- read_rds("./data/G_sila_clean_full_dat.RDS") %>%
  mutate(month_sex = factor(paste(month, sex_M_F),
                            levels = c("April Female", "April Male",
                                       "May Female", "May Male",
                                       "July Female", "July Male")),
         date_chunk_pretty = factor(case_when(month == "April" ~ "Apr 23-25",
                                              month == "May" ~ "May 7-9",
                                              month == "July" ~ "Jul 14"),
                                    levels = c("Apr 23-25", "May 7-9", "Jul 14")),
         gravid_Y_N = factor(gravid_Y_N, levels = c("Not Gravid", "Gravid")),
         percent_water_mass_change = (tmt_mass_change_g/mass_g)*100)
summary(dat)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 10:39:00.00   Min.   :2021-04-23   Length:119         
##  1st Qu.:2021-04-23 15:39:15.00   1st Qu.:2021-04-24   Class :character   
##  Median :2021-04-24 12:22:00.00   Median :2021-04-24   Mode  :character   
##  Mean   :2021-04-28 14:08:18.98   Mean   :2021-05-08                      
##  3rd Qu.:2021-05-07 13:16:45.00   3rd Qu.:2021-05-08                      
##  Max.   :2021-05-08 14:19:00.00   Max.   :2021-07-14                      
##  NA's   :21                                                               
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:119         F-12   :  3   Female:44   Min.   :20.00   Min.   : 77.0  
##  Class :character   M-09   :  3   Male  :75   1st Qu.:31.00   1st Qu.: 93.5  
##  Mode  :character   M-10   :  3               Median :37.00   Median : 99.0  
##                     M-11   :  3               Mean   :36.94   Mean   : 99.0  
##                     M-19   :  3               3rd Qu.:43.00   3rd Qu.:106.5  
##                     M-20   :  3               Max.   :58.80   Max.   :122.0  
##                     (Other):101                                              
##         tmt     tmt_mass_change_g capture_to_msmt   hematocrit_percent
##  Water Tmt:38   Min.   :-3.0000   Min.   :  10.00   Min.   :10.00     
##  Sham Tmt :38   1st Qu.:-0.1500   1st Qu.:  49.25   1st Qu.:30.00     
##  No Tmt   :43   Median : 0.0000   Median : 117.50   Median :34.00     
##                 Mean   :-0.1327   Mean   : 276.18   Mean   :33.78     
##                 3rd Qu.: 0.0000   3rd Qu.: 228.75   3rd Qu.:38.00     
##                 Max.   : 3.0000   Max.   :1665.00   Max.   :58.00     
##                 NA's   :15        NA's   :21        NA's   :3         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :22.00   Min.   : 0.650   Min.   :19.00   Min.   :11.84  
##  1st Qu.:31.00   1st Qu.: 8.486   1st Qu.:28.57   1st Qu.:14.64  
##  Median :33.00   Median :10.381   Median :30.38   Median :16.12  
##  Mean   :32.22   Mean   :10.272   Mean   :29.43   Mean   :20.90  
##  3rd Qu.:34.00   3rd Qu.:12.881   3rd Qu.:31.52   3rd Qu.:23.78  
##  Max.   :36.00   Max.   :21.673   Max.   :33.55   Max.   :40.22  
##                  NA's   :1        NA's   :1       NA's   :1      
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.196   Min.   :1.371   Min.   :309.3           April:67   1 :67  
##  1st Qu.:3.906   1st Qu.:2.946   1st Qu.:348.2           May  :37   3 :37  
##  Median :4.336   Median :3.635   Median :366.5           July :15   12:15  
##  Mean   :4.160   Mean   :3.332   Mean   :367.8                             
##  3rd Qu.:4.626   3rd Qu.:3.960   3rd Qu.:379.9                             
##  Max.   :5.188   Max.   :4.319   Max.   :437.3                             
##  NA's   :1       NA's   :1       NA's   :3                                 
##     week_num      ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   : 1.000   Min.   :2021-04-24   Not Gravid:12   Min.   :3.000  
##  1st Qu.: 1.000   1st Qu.:2021-04-24   Gravid    :16   1st Qu.:3.750  
##  Median : 1.000   Median :2021-04-24   NA's      :91   Median :4.000  
##  Mean   : 3.008   Mean   :2021-04-29                   Mean   :3.812  
##  3rd Qu.: 3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :12.000   Max.   :2021-05-08                   Max.   :5.000  
##                   NA's   :91                           NA's   :103    
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage    
##  Min.   :0.5700          Min.   :1.710                small round:  1  
##  1st Qu.:0.8425          1st Qu.:2.385                large round:  9  
##  Median :0.9950          Median :2.885                soft       :  2  
##  Mean   :1.1131          Mean   :3.148                soft oblong:  1  
##  3rd Qu.:1.4200          3rd Qu.:3.985                firm oblong:  2  
##  Max.   :1.8000          Max.   :4.740                hard oblong:  1  
##  NA's   :103             NA's   :103                  NA's       :103  
##    dev_point         SMI               month_sex  date_chunk_pretty
##  Min.   :1.00   Min.   :22.35   April Female:26   Apr 23-25:67     
##  1st Qu.:2.00   1st Qu.:31.58   April Male  :41   May 7-9  :37     
##  Median :2.00   Median :34.29   May Female  :15   Jul 14   :15     
##  Mean   :2.75   Mean   :34.57   May Male    :22                    
##  3rd Qu.:3.25   3rd Qu.:37.33   July Female : 3                    
##  Max.   :5.00   Max.   :45.63   July Male   :12                    
##  NA's   :103                                                       
##  percent_water_mass_change
##  Min.   :-6.8337          
##  1st Qu.:-0.4487          
##  Median : 0.0000          
##  Mean   :-0.3281          
##  3rd Qu.: 0.0000          
##  Max.   :10.0000          
##  NA's   :15
unique(dat$radio_collar_serial)
##  [1] "245-761"  "252-886"  "245-762"  "252-887"  "245-753"  "245-752" 
##  [7] "245-754"  "229-079"  "245-744"  "229-069"  "252-885"  "245-759" 
## [13] "245-751"  "245-757"  "229-059"  "229-063"  "245-747"  "229-070" 
## [19] "229-065"  "245-741"  "245-748"  "252-882"  "252-881"  "229-078" 
## [25] "245-750"  "229-080"  "245-760"  "245-746"  "245-745"  "245-742" 
## [31] "252-884"  "252-883"  "245-756"  "229-073"  "229-066"  "229-072" 
## [37] "229-060"  "229-074"  "229-077"  "245-762a" "229-070a" NA
dat_F <- dat %>%
  dplyr::filter(sex_M_F == "Female")
dat_M <- dat %>%
  dplyr::filter(sex_M_F == "Male")

# summary statistics
dat %>%
  summarise(mean(CEWL_g_m2h, na.rm = T), sd(CEWL_g_m2h, na.rm = T),
            mean(osmolality_mmol_kg_mean, na.rm = T), 
            sd(osmolality_mmol_kg_mean, na.rm = T),
            mean(hematocrit_percent, na.rm = T), sd(hematocrit_percent, na.rm = T),
            mean(mass_g), sd(mass_g),
            mean(SVL_mm), sd(SVL_mm),
            mean(SMI), sd(SMI))
##   mean(CEWL_g_m2h, na.rm = T) sd(CEWL_g_m2h, na.rm = T)
## 1                    10.27202                  3.924563
##   mean(osmolality_mmol_kg_mean, na.rm = T)
## 1                                 367.8477
##   sd(osmolality_mmol_kg_mean, na.rm = T) mean(hematocrit_percent, na.rm = T)
## 1                               26.73921                            33.77586
##   sd(hematocrit_percent, na.rm = T) mean(mass_g) sd(mass_g) mean(SVL_mm)
## 1                          7.363938     36.94034   8.972626           99
##   sd(SVL_mm) mean(SMI)  sd(SMI)
## 1   9.401587  34.56663 4.475429

Variable Explanation

This is an explanation of each variable in our dataframe:

  • “capture_date_time” = date and time the lizard was captured
  • “capture_date” = date of capture
  • “radio_collar_serial” = serial number of the radio-collar that lizard was fitted with, if applicable; only some lizards were fitted with collars
  • “PIT_tag_ID” = serial number of the PIT tag injected into the lizard at initial capture; all lizards got PIT-tagged
  • “individual_ID” = ID number assigned to lizards
  • “sex_M_F” = lizard sex
  • “mass_g” = lizard body mass
  • “SVL_mm” = lizard snout-vent length
  • “tmt” = treatment given to the lizard: water supplementation, sham, or none (for the lizards that were not radio-collared)
  • “tmt_mass_change_g” = the change in mass for lizards from before to after treatment; only for lizards given water or sham treatments
  • “capture_to_msmt” = time between a lizard being captured and being measured
  • “hematocrit_percent” = percent of blood sample that was red blood cells
  • “cloacal_temp_C” = internal body temperature of the lizard at the time CEWL was taken, taken with a thermocouple in the cloaca
  • “CEWL_g_m2h” = cutaneous evaporative water loss; see wrangling_CEWL for data cleaning
  • “msmt_temp_C” = ambient temperature at the time of CEWL measurement
  • “msmt_RH_percent” = ambient humidity at the time of CEWL measurement
  • “e_s_kPa” = saturation vapor pressure for that ambient temperature
  • “msmt_VPD_kPa” = ambient vapor pressure deficit at the time of CEWL measurement
  • “osmolality_mmol_kg_mean” = plasma osmolality; see wrangling_osml for data cleaning
  • “month” = measurement period as categorical month (April, May, June)
  • “week” = measurement period as categorical week throughout the summer (1, 3, 12)
  • “week_num” = measurement period as numerical week
  • “ultrasound_date” = date that a female was ultrasounded; may differ from capture and measurement dates
  • “gravid_Y_N” = whether or not a female was measured with eggs, if measured (not measured for the non-radio-collared ones)
  • “number_eggs” = if female is gravid, the number of eggs visible with the ultrasound
  • “largest_egg_diameter_cm” = if female is gravid, the diameter of the largest egg in the clutch
  • “largest_egg_circumference_cm” = if female is gravid, the circumference of the largest egg in the clutch
  • “stage” = if female is gravid, approximate stage of development of the eggs
  • “dev_point” = numerical version of development stage
  • “SMI” = body condition of lizards as scaled mass index (calculated in wrangling_general)
  • “month_sex” = variable that I created to make categories for each sex each month
  • “date_chunk_pretty” = nicely formatted measurement period labels

Repeat Subset

I want to make some sub-dataframes for when I want to do stats for only a subset of individuals.

First, get a dataframe of only individuals that had more than one measurement taken, and note how many individuals had each n measurements:

# how many lizards were measured 1, 2, or 3 times?
note_repeats <- dat %>%
  group_by(individual_ID) %>%
  summarise(n = n()) %>%
  group_by(n) %>%
  summarise(n_lizards = n())
note_repeats 
## # A tibble: 3 × 2
##       n n_lizards
##   <int>     <int>
## 1     1        45
## 2     2        28
## 3     3         6
# this is how many individual lizards we measured :)
sum(note_repeats$n_lizards) 
## [1] 79
# get subsetted data
dat_repeats <- dat %>%
  group_by(individual_ID) %>%
  mutate(n = n()) %>%
  dplyr::filter(n>1)

# check these match 'note_repeats'
dat_repeats %>%
  group_by(individual_ID) %>%
  summarise(n = n()) %>%
  group_by(n) %>%
  summarise(n_lizards = n())
## # A tibble: 2 × 2
##       n n_lizards
##   <int>     <int>
## 1     2        28
## 2     3         6

Female Subset

To look at reproductive effort ~ water balance, I will make another subset for only the females we collected gravidity/egg data on:

dat_repro_females <- dat %>%
  dplyr::filter(sex_M_F == "Female") %>%
  dplyr::filter(complete.cases(gravid_Y_N)) %>%
  group_by(individual_ID) %>%
  mutate(n = n(),
         n = factor(n)) 
summary(dat_repro_females)
##  capture_date_time                capture_date        radio_collar_serial
##  Min.   :2021-04-23 11:19:00.0   Min.   :2021-04-23   Length:28          
##  1st Qu.:2021-04-23 15:30:00.0   1st Qu.:2021-04-23   Class :character   
##  Median :2021-04-24 11:11:00.0   Median :2021-04-24   Mode  :character   
##  Mean   :2021-04-28 16:34:38.4   Mean   :2021-04-29                      
##  3rd Qu.:2021-05-08 08:54:00.0   3rd Qu.:2021-05-08                      
##  Max.   :2021-05-08 14:19:00.0   Max.   :2021-05-08                      
##  NA's   :3                                                               
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:28          F-01   : 2    Female:28   Min.   :26.00   Min.   : 85.0  
##  Class :character   F-03   : 2    Male  : 0   1st Qu.:33.52   1st Qu.: 95.0  
##  Mode  :character   F-05   : 2                Median :38.00   Median :100.0  
##                     F-06   : 2                Mean   :39.79   Mean   :101.1  
##                     F-11   : 2                3rd Qu.:44.50   3rd Qu.:108.0  
##                     F-12   : 2                Max.   :58.80   Max.   :114.0  
##                     (Other):16                                               
##         tmt     tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:15   Min.   :-3.0000   Min.   :  25.0   Min.   :18.00     
##  Sham Tmt :13   1st Qu.:-1.0000   1st Qu.:  49.0   1st Qu.:26.75     
##  No Tmt   : 0   Median :-0.1000   Median : 118.0   Median :32.00     
##                 Mean   :-0.4519   Mean   : 238.8   Mean   :32.14     
##                 3rd Qu.: 0.0000   3rd Qu.: 163.0   3rd Qu.:37.00     
##                 Max.   : 1.9000   Max.   :1665.0   Max.   :44.00     
##                 NA's   :1         NA's   :3                          
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :28.00   Min.   : 1.493   Min.   :23.13   Min.   :12.20  
##  1st Qu.:33.00   1st Qu.: 9.171   1st Qu.:30.36   1st Qu.:14.38  
##  Median :34.00   Median :10.239   Median :31.40   Median :15.55  
##  Mean   :33.46   Mean   :11.423   Mean   :30.70   Mean   :16.98  
##  3rd Qu.:35.00   3rd Qu.:13.613   3rd Qu.:31.93   3rd Qu.:19.05  
##  Max.   :36.00   Max.   :21.673   Max.   :33.55   Max.   :28.17  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.830   Min.   :2.033   Min.   :309.3           April:17   1 :17  
##  1st Qu.:4.331   1st Qu.:3.600   1st Qu.:337.8           May  :11   3 :11  
##  Median :4.595   Median :3.863   Median :346.6           July : 0   12: 0  
##  Mean   :4.444   Mean   :3.705   Mean   :352.6                             
##  3rd Qu.:4.736   3rd Qu.:4.050   3rd Qu.:369.2                             
##  Max.   :5.188   Max.   :4.319   Max.   :402.0                             
##                                                                            
##     week_num     ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1.000   Min.   :2021-04-24   Not Gravid:12   Min.   :3.000  
##  1st Qu.:1.000   1st Qu.:2021-04-24   Gravid    :16   1st Qu.:3.750  
##  Median :1.000   Median :2021-04-24                   Median :4.000  
##  Mean   :1.786   Mean   :2021-04-29                   Mean   :3.812  
##  3rd Qu.:3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3.000   Max.   :2021-05-08                   Max.   :5.000  
##                                                       NA's   :12     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.5700          Min.   :1.710                small round: 1  
##  1st Qu.:0.8425          1st Qu.:2.385                large round: 9  
##  Median :0.9950          Median :2.885                soft       : 2  
##  Mean   :1.1131          Mean   :3.148                soft oblong: 1  
##  3rd Qu.:1.4200          3rd Qu.:3.985                firm oblong: 2  
##  Max.   :1.8000          Max.   :4.740                hard oblong: 1  
##  NA's   :12              NA's   :12                   NA's       :12  
##    dev_point         SMI               month_sex  date_chunk_pretty
##  Min.   :1.00   Min.   :27.04   April Female:17   Apr 23-25:17     
##  1st Qu.:2.00   1st Qu.:32.45   April Male  : 0   May 7-9  :11     
##  Median :2.00   Median :34.89   May Female  :11   Jul 14   : 0     
##  Mean   :2.75   Mean   :35.82   May Male    : 0                    
##  3rd Qu.:3.25   3rd Qu.:40.19   July Female : 0                    
##  Max.   :5.00   Max.   :45.63   July Male   : 0                    
##  NA's   :12                                                        
##  percent_water_mass_change n     
##  Min.   :-6.834            2:18  
##  1st Qu.:-3.046            1:10  
##  Median :-0.250                  
##  Mean   :-1.171                  
##  3rd Qu.: 0.000                  
##  Max.   : 7.037                  
##  NA's   :1

It’s good to note that I measured reproductive effort at two time points for 9 females, and only at one time point for 10 females.

Check whether the single-measurement females are all from the first weekend or not:

table(dat_repro_females$week, dat_repro_females$n)
##     
##      2 1
##   1  9 8
##   3  9 2
##   12 0 0
dat_repro_females %>%
  dplyr::filter(n == 1) %>%
  summary() # 4 g, 6 ng
##  capture_date_time              capture_date        radio_collar_serial
##  Min.   :2021-04-23 11:19:00   Min.   :2021-04-23   Length:10          
##  1st Qu.:2021-04-24 09:57:00   1st Qu.:2021-04-24   Class :character   
##  Median :2021-04-24 11:10:00   Median :2021-04-24   Mode  :character   
##  Mean   :2021-04-25 20:16:20   Mean   :2021-04-26                      
##  3rd Qu.:2021-04-24 12:56:00   3rd Qu.:2021-04-24                      
##  Max.   :2021-05-08 14:19:00   Max.   :2021-05-08                      
##  NA's   :1                                                             
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm      
##  Length:10          F-02   :1     Female:10   Min.   :26.00   Min.   : 85.00  
##  Class :character   F-04   :1     Male  : 0   1st Qu.:32.00   1st Qu.: 96.25  
##  Mode  :character   F-07   :1                 Median :37.50   Median :103.50  
##                     F-08   :1                 Mean   :36.30   Mean   :101.80  
##                     F-08a  :1                 3rd Qu.:41.75   3rd Qu.:110.25  
##                     F-09   :1                 Max.   :47.00   Max.   :114.00  
##                     (Other):4                                                 
##         tmt    tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:7   Min.   :-1.00     Min.   :  35.0   Min.   :18.00     
##  Sham Tmt :3   1st Qu.:-1.00     1st Qu.:  49.0   1st Qu.:25.00     
##  No Tmt   :0   Median :-0.50     Median :  91.0   Median :26.00     
##                Mean   :-0.21     Mean   : 261.7   Mean   :29.20     
##                3rd Qu.: 0.00     3rd Qu.: 132.0   3rd Qu.:35.75     
##                Max.   : 1.90     Max.   :1665.0   Max.   :41.00     
##                                  NA's   :1                          
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :28.00   Min.   : 8.457   Min.   :25.03   Min.   :14.12  
##  1st Qu.:32.25   1st Qu.: 9.554   1st Qu.:29.50   1st Qu.:15.05  
##  Median :33.50   Median :10.159   Median :30.96   Median :16.91  
##  Mean   :33.00   Mean   :11.986   Mean   :30.54   Mean   :18.24  
##  3rd Qu.:34.00   3rd Qu.:14.551   3rd Qu.:32.12   3rd Qu.:20.05  
##  Max.   :35.00   Max.   :18.743   Max.   :33.55   Max.   :26.03  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month   week  
##  Min.   :3.172   Min.   :2.346   Min.   :309.3           April:8   1 :8  
##  1st Qu.:4.124   1st Qu.:3.341   1st Qu.:330.6           May  :2   3 :2  
##  Median :4.482   Median :3.776   Median :357.3           July :0   12:0  
##  Mean   :4.410   Mean   :3.623   Mean   :355.6                           
##  3rd Qu.:4.786   3rd Qu.:4.014   3rd Qu.:378.2                           
##  Max.   :5.188   Max.   :4.319   Max.   :402.0                           
##                                                                          
##     week_num   ultrasound_date           gravid_Y_N  number_eggs  
##  Min.   :1.0   Min.   :2021-04-24   Not Gravid:6    Min.   :3.00  
##  1st Qu.:1.0   1st Qu.:2021-04-24   Gravid    :4    1st Qu.:3.75  
##  Median :1.0   Median :2021-04-24                   Median :4.00  
##  Mean   :1.4   Mean   :2021-04-26                   Mean   :3.75  
##  3rd Qu.:1.0   3rd Qu.:2021-04-24                   3rd Qu.:4.00  
##  Max.   :3.0   Max.   :2021-05-08                   Max.   :4.00  
##                                                     NA's   :6     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage  
##  Min.   :0.570           Min.   :1.710                small round:0  
##  1st Qu.:0.795           1st Qu.:2.265                large round:3  
##  Median :0.935           Median :2.595                soft       :0  
##  Mean   :0.935           Mean   :2.717                soft oblong:0  
##  3rd Qu.:1.075           3rd Qu.:3.047                firm oblong:1  
##  Max.   :1.300           Max.   :3.970                hard oblong:0  
##  NA's   :6               NA's   :6                    NA's       :6  
##    dev_point         SMI               month_sex date_chunk_pretty
##  Min.   :2.00   Min.   :27.04   April Female:8   Apr 23-25:8      
##  1st Qu.:2.00   1st Qu.:30.81   April Male  :0   May 7-9  :2      
##  Median :2.00   Median :32.93   May Female  :2   Jul 14   :0      
##  Mean   :2.75   Mean   :32.21   May Male    :0                    
##  3rd Qu.:2.75   3rd Qu.:33.44   July Female :0                    
##  Max.   :5.00   Max.   :37.68   July Male   :0                    
##  NA's   :6                                                        
##  percent_water_mass_change n     
##  Min.   :-3.8462           2: 0  
##  1st Qu.:-2.6849           1:10  
##  Median :-1.1628                 
##  Mean   :-0.5466                 
##  3rd Qu.: 0.0000                 
##  Max.   : 7.0370                 
## 
dat_repro_females %>%
  dplyr::filter(n == 2) %>%
  dplyr::filter(week == 3) %>%
  summary() # 8 g, 1 ng
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-05-08 08:40:00.00   Min.   :2021-05-08   Length:9           
##  1st Qu.:2021-05-08 09:09:00.00   1st Qu.:2021-05-08   Class :character   
##  Median :2021-05-08 09:58:00.00   Median :2021-05-08   Mode  :character   
##  Mean   :2021-05-08 10:06:42.85   Mean   :2021-05-08                      
##  3rd Qu.:2021-05-08 10:42:00.00   3rd Qu.:2021-05-08                      
##  Max.   :2021-05-08 12:27:00.00   Max.   :2021-05-08                      
##  NA's   :2                                                                
##   PIT_tag_ID        individual_ID   sex_M_F      mass_g          SVL_mm      
##  Length:9           F-01   :1     Female:9   Min.   :33.70   Min.   : 90.00  
##  Class :character   F-03   :1     Male  :0   1st Qu.:35.60   1st Qu.: 95.00  
##  Mode  :character   F-05   :1                Median :43.90   Median : 98.00  
##                     F-06   :1                Mean   :44.02   Mean   : 99.56  
##                     F-11   :1                3rd Qu.:52.70   3rd Qu.:105.00  
##                     F-12   :1                Max.   :58.80   Max.   :110.00  
##                     (Other):3                                                
##         tmt    tmt_mass_change_g capture_to_msmt hematocrit_percent
##  Water Tmt:4   Min.   :-3.0      Min.   : 59.0   Min.   :31.00     
##  Sham Tmt :5   1st Qu.:-1.5      1st Qu.:158.5   1st Qu.:32.00     
##  No Tmt   :0   Median :-1.0      Median :244.0   Median :37.00     
##                Mean   :-1.0      Mean   :201.7   Mean   :37.33     
##                3rd Qu.: 0.0      3rd Qu.:259.5   3rd Qu.:43.00     
##                Max.   : 1.0      Max.   :273.0   Max.   :44.00     
##                NA's   :1         NA's   :2                         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :33.00   Min.   : 8.318   Min.   :30.40   Min.   :13.12  
##  1st Qu.:34.00   1st Qu.: 9.318   1st Qu.:30.86   1st Qu.:14.44  
##  Median :35.00   Median :10.286   Median :31.45   Median :14.70  
##  Mean   :34.56   Mean   :10.915   Mean   :31.38   Mean   :14.78  
##  3rd Qu.:35.00   3rd Qu.:12.213   3rd Qu.:31.90   3rd Qu.:15.48  
##  Max.   :36.00   Max.   :14.378   Max.   :32.04   Max.   :15.84  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month   week  
##  Min.   :4.340   Min.   :3.685   Min.   :337.3           April:0   1 :0  
##  1st Qu.:4.456   1st Qu.:3.813   1st Qu.:340.5           May  :9   3 :9  
##  Median :4.608   Median :3.939   Median :361.0           July :0   12:0  
##  Mean   :4.591   Mean   :3.913   Mean   :357.0                           
##  3rd Qu.:4.728   3rd Qu.:4.028   3rd Qu.:369.0                           
##  Max.   :4.765   Max.   :4.112   Max.   :379.5                           
##                                                                          
##     week_num ultrasound_date           gravid_Y_N  number_eggs  
##  Min.   :3   Min.   :2021-05-08   Not Gravid:1    Min.   :3.00  
##  1st Qu.:3   1st Qu.:2021-05-08   Gravid    :8    1st Qu.:3.00  
##  Median :3   Median :2021-05-08                   Median :4.00  
##  Mean   :3   Mean   :2021-05-08                   Mean   :3.75  
##  3rd Qu.:3   3rd Qu.:2021-05-08                   3rd Qu.:4.00  
##  Max.   :3   Max.   :2021-05-08                   Max.   :5.00  
##                                                   NA's   :1     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage  
##  Min.   :0.640           Min.   :1.820                small round:1  
##  1st Qu.:1.135           1st Qu.:3.333                large round:2  
##  Median :1.440           Median :4.000                soft       :2  
##  Mean   :1.335           Mean   :3.712                soft oblong:1  
##  3rd Qu.:1.560           3rd Qu.:4.228                firm oblong:1  
##  Max.   :1.800           Max.   :4.740                hard oblong:1  
##  NA's   :1               NA's   :1                    NA's       :1  
##    dev_point          SMI               month_sex date_chunk_pretty
##  Min.   :1.000   Min.   :34.49   April Female:0   Apr 23-25:0      
##  1st Qu.:2.000   1st Qu.:36.44   April Male  :0   May 7-9  :9      
##  Median :3.000   Median :41.85   May Female  :9   Jul 14   :0      
##  Mean   :3.125   Mean   :40.73   May Male    :0                    
##  3rd Qu.:4.250   3rd Qu.:43.87   July Female :0                    
##  Max.   :5.000   Max.   :45.63   July Male   :0                    
##  NA's   :1                                                         
##  percent_water_mass_change n    
##  Min.   :-6.834            2:9  
##  1st Qu.:-3.649            1:0  
##  Median :-2.468                 
##  Mean   :-2.203                 
##  3rd Qu.: 0.000                 
##  Max.   : 2.809                 
##  NA's   :1
table(dat_repro_females$gravid_Y_N, dat_repro_females$n)
##             
##               2  1
##   Not Gravid  6  6
##   Gravid     12  4

For the females who were only measured once, 8 of them are from the first weekend and 2 of them are from the third weekend.

Also save a subset of the females that were measured twice:

dat_repeat_females <- dat_repro_females %>%
  dplyr::filter(n == 2)

Does SVL change over time?

I want to know whether SVL changed over time for the individuals who had more than one measurement taken.

# model
SVL_change <- lmerTest::lmer(data = dat_repeats, 
                   SVL_mm ~ week_num + (1|individual_ID))
summary(SVL_change)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SVL_mm ~ week_num + (1 | individual_ID)
##    Data: dat_repeats
## 
## REML criterion at convergence: 495.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.46209 -0.38905 -0.04905  0.51484  2.46481 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_ID (Intercept) 55.65    7.460   
##  Residual                  19.96    4.467   
## Number of obs: 74, groups:  individual_ID, 34
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 100.5373     1.4817  42.0626  67.852   <2e-16 ***
## week_num      0.1016     0.1442  42.5756   0.705    0.485    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## week_num -0.359

There is no meaningful change in SVL for the individuals we were able to measure, so I could use either mass or body condition to estimate change in body size over time.

Check CEWL Covariates

How much is CEWL impacted by the body temperature, ambient temperature, and ambient VPD at the time of CEWL measurement? I will investigate this using all measurements.

It’s possible that I may find a relationship between ambient temp and CEWL but not body temp and CEWL because the body temp measurements we took were unfortunately only to the nearest degree (vs tenth of a degree).

Plot

First, make a simple plot of each to look for relationships visually.

ggplot(dat) +
  aes(x = cloacal_temp_C,
      y = CEWL_g_m2h,
      color = week
      #color = sex_M_F
        ) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw()

ggplot(dat) +
  aes(x = msmt_temp_C,
      y = CEWL_g_m2h,
      color = week
      #color = sex_M_F
        ) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw()

ggplot(dat) +
  aes(x = msmt_VPD_kPa,
      y = CEWL_g_m2h,
      color = week
      #color = sex_M_F
        ) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw()

Across dates, CEWL had no relationship with body temp, ambient temp, or ambient VPD at the time of capture. However, when separated by date, there were relationships. CEWL had a positive relationship with body temp, ambient temp, and ambient VPD, but only for week 3 measurements. There were no relationships when separated by treatment group.

Models

Now, check linear models of the effect of each of body temp, ambient temp, and ambient VPD with the interaction of date on CEWL.

CEWL_cov_check_clotemp <- lm(data = dat,
           CEWL_g_m2h ~ cloacal_temp_C*week)
summary(CEWL_cov_check_clotemp)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ cloacal_temp_C * week, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1993  -2.0150   0.1478   2.0471  10.8241 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            11.45021    4.57266   2.504 0.013718 *  
## cloacal_temp_C         -0.01821    0.14195  -0.128 0.898154    
## week3                 -63.35695   17.24290  -3.674 0.000367 ***
## week12                -20.29849   29.16728  -0.696 0.487912    
## cloacal_temp_C:week3    1.82229    0.51531   3.536 0.000591 ***
## cloacal_temp_C:week12   0.71647    0.98516   0.727 0.468586    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.613 on 112 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1888, Adjusted R-squared:  0.1526 
## F-statistic: 5.213 on 5 and 112 DF,  p-value: 0.0002434
CEWL_cov_check_msmttemp <- lm(data = dat,
           CEWL_g_m2h ~ msmt_temp_C*week)
summary(CEWL_cov_check_msmttemp)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ msmt_temp_C * week, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2517  -1.5583  -0.3139   2.0494  10.7491 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         10.25244    3.27104   3.134   0.0022 ** 
## msmt_temp_C          0.02113    0.11155   0.189   0.8501    
## week3              -60.68682   14.33616  -4.233 4.74e-05 ***
## week12              -3.76401   34.56181  -0.109   0.9135    
## msmt_temp_C:week3    1.90326    0.46806   4.066 8.91e-05 ***
## msmt_temp_C:week12   0.16788    1.23504   0.136   0.8921    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.554 on 112 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2149, Adjusted R-squared:  0.1799 
## F-statistic: 6.133 on 5 and 112 DF,  p-value: 4.648e-05
CEWL_cov_check_VPD <- lm(data = dat,
           CEWL_g_m2h ~ msmt_VPD_kPa*week)
summary(CEWL_cov_check_VPD)
## 
## Call:
## lm(formula = CEWL_g_m2h ~ msmt_VPD_kPa * week, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1892  -1.7946  -0.3724   2.0215  10.8386 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          11.13212    1.80429   6.170 1.12e-08 ***
## msmt_VPD_kPa         -0.07911    0.52156  -0.152    0.880    
## week3               -25.66036    5.31654  -4.827 4.43e-06 ***
## week12                0.33315   13.99194   0.024    0.981    
## msmt_VPD_kPa:week3    6.37242    1.44871   4.399 2.49e-05 ***
## msmt_VPD_kPa:week12   0.20939    5.87420   0.036    0.972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.504 on 112 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2368, Adjusted R-squared:  0.2027 
## F-statistic:  6.95 on 5 and 112 DF,  p-value: 1.093e-05
# F stats
car::Anova(CEWL_cov_check_clotemp, type = "II")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##                      Sum Sq  Df F value   Pr(>F)   
## cloacal_temp_C        12.30   1  0.9422 0.333793   
## week                 169.15   2  6.4797 0.002173 **
## cloacal_temp_C:week  167.73   2  6.4253 0.002282 **
## Residuals           1461.83 112                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(CEWL_cov_check_msmttemp, type = "II")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##                   Sum Sq  Df F value    Pr(>F)    
## msmt_temp_C        18.24   1  1.4441 0.2320199    
## week              178.22   2  7.0546 0.0013014 ** 
## msmt_temp_C:week  208.88   2  8.2684 0.0004474 ***
## Residuals        1414.73 112                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(CEWL_cov_check_VPD, type = "II")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##                    Sum Sq  Df F value    Pr(>F)    
## msmt_VPD_kPa        28.80   1  2.3452 0.1284889    
## week               188.86   2  7.6898 0.0007423 ***
## msmt_VPD_kPa:week  237.73   2  9.6798 0.0001325 ***
## Residuals         1375.33 112                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CEWL had a positive relationship with each of body temp, ambient temp, and ambient VPD, but only for the measurements taken in May (week 3).

Variance

Let’s also check the distribution/range of these covariates, because based on the figures, week 3 had a pretty small range of temperatures and VPDs.

dat %>%
  group_by(week) %>%
  summarise(ctemp_range = max(cloacal_temp_C, na.rm = T) - 
              min(cloacal_temp_C, na.rm = T),
            mtemp_range = max(msmt_temp_C, na.rm = T) - 
              min(msmt_temp_C, na.rm = T),
            mVPD_range = max(msmt_VPD_kPa, na.rm = T) - 
              min(msmt_VPD_kPa, na.rm = T))
## # A tibble: 3 × 4
##   week  ctemp_range mtemp_range mVPD_range
##   <fct>       <dbl>       <dbl>      <dbl>
## 1 1              14       14.5       2.95 
## 2 3               5        4.55      1.28 
## 3 12              3        2.22      0.495

Well… the ranges for week 3 are not small enough to be meaningless. So, we will have to keep in mind that temperature and VPD should be retained in CEWL models as a covariate.

Water Balance Metrics

What is the interrelation among water balance metrics CEWL, osmolality, hematocrit, and mass/SMI. I’ll use all measurements for this.

Models

Check a simple linear model between each set of paired measurements. I started with full models that included interactions with month and sex for each relationship, then dropped if the interaction was not significant. So, the models that still have interaction effects have them because they are important, and the models that do not have interaction effects had those effects dropped because they were not important.

# CEWL ~
car::Anova(lm(data = dat,
           CEWL_g_m2h ~ osmolality_mmol_kg_mean*month), type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##                                Sum Sq  Df F value  Pr(>F)  
## osmolality_mmol_kg_mean         17.59   1  1.3402 0.24950  
## month                          125.33   2  4.7747 0.01027 *
## osmolality_mmol_kg_mean:month   72.58   2  2.7652 0.06734 .
## Residuals                     1443.69 110                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(data = dat,
           CEWL_g_m2h ~ osmolality_mmol_kg_mean*month))
## 
## Call:
## lm(formula = CEWL_g_m2h ~ osmolality_mmol_kg_mean * month, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3779 -2.0003  0.0237  1.9817 10.8246 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                        7.266624   6.658883   1.091   0.2775  
## osmolality_mmol_kg_mean            0.010368   0.018227   0.569   0.5706  
## monthMay                          21.295848  10.315446   2.064   0.0413 *
## monthJuly                         -2.159012  23.011966  -0.094   0.9254  
## osmolality_mmol_kg_mean:monthMay  -0.062972   0.027584  -2.283   0.0244 *
## osmolality_mmol_kg_mean:monthJuly  0.008481   0.064845   0.131   0.8962  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.623 on 110 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.1548, Adjusted R-squared:  0.1164 
## F-statistic:  4.03 on 5 and 110 DF,  p-value: 0.00214
car::Anova(lm(data = dat,
           CEWL_g_m2h ~ hematocrit_percent), type = "II")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##                    Sum Sq  Df F value   Pr(>F)   
## hematocrit_percent  136.0   1  9.8621 0.002149 **
## Residuals          1572.1 114                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(lm(data = dat,
           CEWL_g_m2h ~ mass_g), type = "2")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##            Sum Sq  Df F value Pr(>F)
## mass_g      22.14   1  1.4432 0.2321
## Residuals 1779.91 116
car::Anova(lm(data = dat,
           CEWL_g_m2h ~ SMI), type = "II")
## Anova Table (Type II tests)
## 
## Response: CEWL_g_m2h
##            Sum Sq  Df F value  Pr(>F)  
## SMI         58.54   1  3.8949 0.05081 .
## Residuals 1743.51 116                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# osmolality ~
car::Anova(lm(data = dat,
           osmolality_mmol_kg_mean ~ hematocrit_percent), type = "II")
## Anova Table (Type II tests)
## 
## Response: osmolality_mmol_kg_mean
##                    Sum Sq  Df F value    Pr(>F)    
## hematocrit_percent   9602   1  15.072 0.0001738 ***
## Residuals           72622 114                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(lm(data = dat,
           osmolality_mmol_kg_mean ~ mass_g), type = "II")
## Anova Table (Type II tests)
## 
## Response: osmolality_mmol_kg_mean
##           Sum Sq  Df F value Pr(>F)
## mass_g         0   1   3e-04 0.9869
## Residuals  82223 114
car::Anova(lm(data = dat,
           osmolality_mmol_kg_mean ~ SMI), type = "II")
## Anova Table (Type II tests)
## 
## Response: osmolality_mmol_kg_mean
##           Sum Sq  Df F value  Pr(>F)  
## SMI         3656   1  5.3045 0.02308 *
## Residuals  78568 114                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# hct
car::Anova(lm(data = dat,
           hematocrit_percent ~ SMI*sex_M_F), type = "II")
## Anova Table (Type II tests)
## 
## Response: hematocrit_percent
##             Sum Sq  Df F value   Pr(>F)   
## SMI          397.6   1  8.8942 0.003511 **
## sex_M_F      468.8   1 10.4855 0.001583 **
## SMI:sex_M_F  356.3   1  7.9694 0.005632 **
## Residuals   5007.3 112                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notes on these relationships:

  • CEWL ~ osml had significant interaction of osml* month (osml* mo* sex was not sig)
  • CEWL ~ hct was sig but had no extra interactions (plot direct relationship only)
  • CEWL ~ SMI was sig (mass was not; & SMI did not interact with mo or sex)
  • osml & hct corr, but with no interactions of mo or sex
  • osml also corr with SMI but not mass (& no interactions with mo or sex)
  • ^all of the above based on anova SS, so assessing whether these things explain variation in response (CEWL/osml)

For MS:

  • Plasma osmolality explained a significant amount of the variation in CEWL (SLR: F(1,110) = 5.07, p = 0.03), with a marginally nonsignificant interaction with measurement time period (SLR: F(2,110) = 2.77, p = 0.07).

  • CEWL could also be explained by hematocrit (SLR: F(1,114) = 9.86, p = 0.002), or body condition (SLR: F(1,116) = 3.89, p = 0.05), but not lizard mass (SLR: F(1,116) = 1.44, p = 0.2).

  • Variation in plasma osmolality could be explained by hematocrit (SLR: F(1,114) = 15.07, p < 0.001) or body condition (SLR: F(1,114) = 5.30, p = 0.02), but not lizard mass (SLR: F(1,114) < 0.001, p = 1).

  • Hematocrit and body condition were also significantly related (SLR: F(1,112) = 9.03, p = 0.003), with an interaction of lizard sex (SLR: F(1,112) = 7.97, p = 0.005).

Plots

CEWL ~ Osmolality

ggplot(dat) + 
  aes(x = osmolality_mmol_kg_mean,
      y = CEWL_g_m2h,
      color = month,
      shape = month
      ) +
  geom_point(alpha = 0.5) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              alpha = 0.8) + 
  scale_color_manual(values = n3, name = NULL
                     #name = "Measurement\nPeriod"
                     ) + 
  scale_shape_manual(values = c(15, 16, 17), name = NULL) +
  labs(x = bquote('Plasma Osmolality (mmol '*kg^-1*')'),
       y = bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
  scale_x_continuous(limits = c(300, 450),
                     breaks = c(300, 350, 400, 450),
                     labels = c(300, 350, 400, 450)) +
  savs_ggplot_theme -> CEWL_osml_fig
CEWL_osml_fig

# arrange to get legend centered
ggarrange(CEWL_osml_fig,
          ncol = 1, nrow = 1,
          common.legend = TRUE,
          legend = "bottom"
          ) -> CEWL_osml_fig_arrange

# export figure
ggsave(filename = "CEWL_osml_fig.pdf",
       plot = CEWL_osml_fig_arrange,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 80)

CEWL ~ SMI

ggplot(dat) + 
  aes(x = SMI,
      y = CEWL_g_m2h,
      #color = tmt
      #color = month_sex
      ) +
  geom_point(size = 1, 
             alpha = 0.4) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1,
              color = n1) + 
  xlab("Body Condition (g')") + 
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  savs_ggplot_theme -> CEWL_SMI_fig
CEWL_SMI_fig

CEWL ~ Hematocrit

ggplot(dat) + 
  aes(x = hematocrit_percent,
      y = CEWL_g_m2h,
      #color = tmt
      #color = month_sex
      ) +
  geom_point(size = 1, 
             alpha = 0.4) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1,
              color = n1) + 
  xlab("Hematocrit (%)") + 
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  #ylab("") +
  #xlim(300, 450) +
  #annotate(geom = "text", x = 430, y = 24, label = bquote(''*R^2*'=0.033'), 
   #        size = 7) +
  #ylim(0, 25) +
  savs_ggplot_theme -> CEWL_hct_fig
CEWL_hct_fig

Osmol ~ SMI

ggplot(dat) + 
  aes(x = SMI,
      y = osmolality_mmol_kg_mean,
      #color = tmt
      #color = month_sex
      ) +
  geom_point(size = 1, 
             alpha = 0.4) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1,
              color = n1) + 
  xlab("Body Condition (g')") + 
  ylab("Osmolality (mmol/kg)") + 
  #ylab("") +
  #xlim(300, 450) +
  #annotate(geom = "text", x = 430, y = 24, label = bquote(''*R^2*'=0.033'), 
   #        size = 7) +
  #ylim(0, 25) +
  savs_ggplot_theme -> osml_SMI_fig
osml_SMI_fig

Osmol ~ Hematocrit

ggplot(dat) + 
  aes(x = hematocrit_percent,
      y = osmolality_mmol_kg_mean,
      #color = tmt
      #color = month_sex
      ) +
  geom_point(size = 1, 
             alpha = 0.4) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1,
              color = n1) + 
  xlab("Hematocrit (%)") + 
  ylab("Osmolality (mmol/kg)") + 
  #ylab("") +
  #xlim(300, 450) +
  #annotate(geom = "text", x = 430, y = 24, label = bquote(''*R^2*'=0.033'), 
   #        size = 7) +
  #ylim(0, 25) +
  scale_color_manual(values = n1) + 
  savs_ggplot_theme -> osml_hct_fig
osml_hct_fig

Hct ~ SMI

ggplot(dat) + 
  aes(x = SMI,
      y = hematocrit_percent,
      color = sex_M_F,
      #color = month_sex
      ) +
  geom_point(size = 1, 
             alpha = 0.4) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1.6, 
              alpha = 1) + 
  xlab("Body Condition (g')") + 
  ylab("Hematocrit (%)") + 
  #ylab("") +
  #xlim(300, 450) +
  #annotate(geom = "text", x = 430, y = 24, label = bquote(''*R^2*'=0.033'), 
   #        size = 7) +
  #ylim(0, 25) +
  scale_color_manual(values = n2) +
  savs_ggplot_theme -> hct_SMI_fig
hct_SMI_fig

Supplemental Hydration Treatment Effects

We tried to give the lizards water to drink. Did they actually drink any? And did it affect measurable change in their physiology?

Calculate Deltas

First, calculate the change in CEWL, osmolality, and mass for each treatment individual that was recaptured between April and May.

# get the intitial values from April
starting_values <- dat_repeats %>%
  # look only at lizards in hydration or sham tmt
  dplyr::filter(tmt != "No Tmt") %>%
  # get the values for April measurements only 
  dplyr::filter(month == "April") %>%
  dplyr::select(individual_ID,
                CEWL_April = CEWL_g_m2h,
                osml_April = osmolality_mmol_kg_mean,
                mass_April = mass_g)

# calculate deltas
dat_tmt_change <- dat_repeats %>%
  # look only at lizards in hydration or sham tmt
  dplyr::filter(tmt != "No Tmt") %>%
  # add starting values
  left_join(starting_values, by = 'individual_ID') %>%
  # make sure we have complete data
  dplyr::filter(complete.cases(mass_April, CEWL_April)) %>%
  # calculate difference for each msmt versus in April
  mutate(mass_change_from_April = mass_g - mass_April,
         CEWL_change_from_April = CEWL_g_m2h - CEWL_April,
         osml_change_from_April = osmolality_mmol_kg_mean - osml_April)

# subset of data to look at change April->May only
dat_tmt_change_May_only <- dat_tmt_change %>%
  dplyr::filter(month == "May")

Tmt Effect Check

What was mass change during treatment, when we either gave lizards water to drink or a sham treatment?

This quantifies the difference in mass gain/loss when we gave them water to drink versus sham treatment.

tmt_effect_mass <- dat %>% 
  dplyr::filter(tmt != "No Tmt") %>% 
  dplyr::filter(complete.cases(tmt_mass_change_g))
summary(tmt_effect_mass)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 10:39:00.00   Min.   :2021-04-23   Length:61          
##  1st Qu.:2021-04-23 14:54:00.00   1st Qu.:2021-04-23   Class :character   
##  Median :2021-04-24 10:49:00.00   Median :2021-04-24   Mode  :character   
##  Mean   :2021-04-28 09:55:24.21   Mean   :2021-04-28                      
##  3rd Qu.:2021-05-07 11:31:00.00   3rd Qu.:2021-05-07                      
##  Max.   :2021-05-08 14:19:00.00   Max.   :2021-05-08                      
##  NA's   :4                                                                
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm   
##  Length:61          F-01   : 2    Female:29   Min.   :26.00   Min.   : 85  
##  Class :character   F-05   : 2    Male  :32   1st Qu.:33.00   1st Qu.: 95  
##  Mode  :character   F-06   : 2                Median :38.00   Median : 99  
##                     F-11   : 2                Mean   :39.38   Mean   :101  
##                     F-12   : 2                3rd Qu.:46.00   3rd Qu.:108  
##                     F-13   : 2                Max.   :56.00   Max.   :122  
##                     (Other):49                                             
##         tmt     tmt_mass_change_g capture_to_msmt  hematocrit_percent
##  Water Tmt:29   Min.   :-3.0000   Min.   :  11.0   Min.   :17.00     
##  Sham Tmt :32   1st Qu.:-1.0000   1st Qu.:  62.0   1st Qu.:29.75     
##  No Tmt   : 0   Median : 0.0000   Median : 150.0   Median :36.00     
##                 Mean   :-0.2262   Mean   : 403.5   Mean   :34.30     
##                 3rd Qu.: 0.0000   3rd Qu.: 266.0   3rd Qu.:39.00     
##                 Max.   : 3.0000   Max.   :1665.0   Max.   :58.00     
##                                   NA's   :4        NA's   :1         
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :24.00   Min.   : 0.650   Min.   :19.07   Min.   :12.20  
##  1st Qu.:32.00   1st Qu.: 6.490   1st Qu.:28.57   1st Qu.:14.53  
##  Median :33.00   Median : 9.318   Median :30.50   Median :18.73  
##  Mean   :32.44   Mean   : 9.410   Mean   :29.43   Mean   :19.45  
##  3rd Qu.:34.00   3rd Qu.:10.827   3rd Qu.:31.86   3rd Qu.:22.46  
##  Max.   :36.00   Max.   :21.673   Max.   :33.55   Max.   :35.83  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month    week   
##  Min.   :2.205   Min.   :1.415   Min.   :309.3           April:39   1 :39  
##  1st Qu.:3.907   1st Qu.:3.014   1st Qu.:346.1           May  :22   3 :22  
##  Median :4.365   Median :3.685   Median :364.2           July : 0   12: 0  
##  Mean   :4.168   Mean   :3.396   Mean   :366.5                             
##  3rd Qu.:4.717   3rd Qu.:4.028   3rd Qu.:381.1                             
##  Max.   :5.188   Max.   :4.319   Max.   :436.5                             
##                                  NA's   :1                                 
##     week_num     ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1.000   Min.   :2021-04-24   Not Gravid:12   Min.   :3.000  
##  1st Qu.:1.000   1st Qu.:2021-04-24   Gravid    :15   1st Qu.:3.500  
##  Median :1.000   Median :2021-04-24   NA's      :34   Median :4.000  
##  Mean   :1.721   Mean   :2021-04-29                   Mean   :3.733  
##  3rd Qu.:3.000   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3.000   Max.   :2021-05-08                   Max.   :4.000  
##                  NA's   :34                           NA's   :46     
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage   
##  Min.   :0.570           Min.   :1.710                small round: 1  
##  1st Qu.:0.815           1st Qu.:2.320                large round: 9  
##  Median :0.990           Median :2.880                soft       : 2  
##  Mean   :1.073           Mean   :3.045                soft oblong: 0  
##  3rd Qu.:1.350           3rd Qu.:3.970                firm oblong: 2  
##  Max.   :1.800           Max.   :4.740                hard oblong: 1  
##  NA's   :46              NA's   :46                   NA's       :46  
##    dev_point          SMI               month_sex  date_chunk_pretty
##  Min.   :1.000   Min.   :27.04   April Female:19   Apr 23-25:39     
##  1st Qu.:2.000   1st Qu.:32.99   April Male  :20   May 7-9  :22     
##  Median :2.000   Median :34.76   May Female  :10   Jul 14   : 0     
##  Mean   :2.667   Mean   :35.40   May Male    :12                    
##  3rd Qu.:3.000   3rd Qu.:38.39   July Female : 0                    
##  Max.   :5.000   Max.   :45.63   July Male   : 0                    
##  NA's   :46                                                         
##  percent_water_mass_change
##  Min.   :-6.8337          
##  1st Qu.:-2.7027          
##  Median : 0.0000          
##  Mean   :-0.5594          
##  3rd Qu.: 0.0000          
##  Max.   :10.0000          
## 
# change in grams (actual amount of water "drank")
tmt_change <- lm(data = tmt_effect_mass, tmt_mass_change_g ~ tmt)
summary(tmt_change)
## 
## Call:
## lm(formula = tmt_mass_change_g ~ tmt, data = tmt_effect_mass)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2344 -0.3690 -0.2344  0.6656  3.7656 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.3690     0.2090   1.765 0.082666 .  
## tmtSham Tmt  -1.1346     0.2886  -3.932 0.000224 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.125 on 59 degrees of freedom
## Multiple R-squared:  0.2076, Adjusted R-squared:  0.1942 
## F-statistic: 15.46 on 1 and 59 DF,  p-value: 0.0002241
emmeans(tmt_change, "tmt")
##  tmt       emmean    SE df lower.CL upper.CL
##  Water Tmt  0.369 0.209 59  -0.0492    0.787
##  Sham Tmt  -0.766 0.199 59  -1.1637   -0.368
## 
## Confidence level used: 0.95
# change in percent body mass (relative amount)
percent_change <- lm(data = tmt_effect_mass, percent_water_mass_change ~ tmt)
summary(percent_change)
## 
## Call:
## lm(formula = percent_water_mass_change ~ tmt, data = tmt_effect_mass)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4465 -1.4239 -0.8996  1.8787 11.9094 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.9303     0.5617   1.656 0.102997    
## tmtSham Tmt  -2.8397     0.7756  -3.661 0.000538 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.025 on 59 degrees of freedom
## Multiple R-squared:  0.1852, Adjusted R-squared:  0.1713 
## F-statistic: 13.41 on 1 and 59 DF,  p-value: 0.0005376
emmeans(percent_change, "tmt")
##  tmt       emmean    SE df lower.CL upper.CL
##  Water Tmt   0.93 0.562 59   -0.194    2.054
##  Sham Tmt   -1.91 0.535 59   -2.979   -0.839
## 
## Confidence level used: 0.95

So, it did result in a significant change in mass within the same-day water supplementation tmt (vs sham), but only ~1.4g difference, on average.

Change Over Time

Did hydric physiology change throughout these lizards’ active season? Based on hydration treatment? Sex?

Plot Change ~ Tmt

Mass

ggplot(dat_tmt_change) +
  aes(x = month,
      y = mass_change_from_April,
      group = individual_ID,
      color = tmt
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'Lizard Mass (g)')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  savs_ggplot_theme -> tmt_delta_mass
tmt_delta_mass

Osmolality

ggplot(dat_tmt_change) +
  aes(x = month,
      y = osml_change_from_April,
      group = individual_ID,
      color = tmt
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'Osmolality (mmol '*kg^-1*')')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  savs_ggplot_theme -> tmt_delta_osml
tmt_delta_osml

CEWL

ggplot(dat_tmt_change) +
  aes(x = month,
      y = CEWL_change_from_April,
      group = individual_ID,
      color = tmt
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'CEWL (g '*m^-2*' '*h^-1*')')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  savs_ggplot_theme -> tmt_delta_CEWL
tmt_delta_CEWL

It does not look like treatment led to any differences in the amount or direction of change in lizard mass, CEWL, or osmolality. Changes were also not consistent by sex.

Export Multi-Panel Fig

Save plots:

ggarrange(tmt_delta_CEWL,
          tmt_delta_osml, 
          tmt_delta_mass,
          ncol = 1, nrow = 3,
          labels = c("A", "B", "C"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> tmt_delta_multi_fig
tmt_delta_multi_fig

Tmt Marginal Means

I want to know whether the average values for the lizards we measured changes over time based on supplemental hydration treatment. To test this, I can run a linear model, then use emmeans to get the marginal means and confidence intervals for each time point.

Do this only for lizards with repeat measurements, because those are the ones that received tmt.

Mass

tmt_diffs_mass <- lmerTest::lmer(data = dat_repeats, mass_g ~ month*tmt +
                                        (1|individual_ID))
anova(tmt_diffs_mass, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##            Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## month     274.074 137.037     2 37.802  7.4633 0.001854 **
## tmt       224.282 112.141     2 32.029  6.1078 0.005661 **
## month:tmt  37.189  12.396     3 37.177  0.6751 0.572776   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
(emmeans(tmt_diffs_mass, pairwise ~ month | tmt)$emmean)
## tmt = Water Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   41.0 2.03 45.7     37.0     45.1
##  May     46.1 2.26 56.7     41.6     50.6
##  July    37.6 2.34 59.7     32.9     42.3
## 
## tmt = Sham Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   37.9 1.92 43.5     34.1     41.8
##  May     39.7 1.99 47.8     35.7     43.7
##  July    34.9 2.53 65.4     29.9     40.0
## 
## tmt = No Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   29.8 3.33 43.5     23.1     36.5
##  May     31.8 3.33 43.5     25.1     38.5
##  July  nonEst   NA   NA       NA       NA
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
# stat diffs by month
(emmeans(tmt_diffs_mass, pairwise ~ month | tmt)$contrasts)
## tmt = Water Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -5.05 1.99 37.6  -2.542  0.0396
##  April - July     3.45 2.09 38.3   1.652  0.2368
##  May - July       8.50 2.44 41.5   3.480  0.0033
## 
## tmt = Sham Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -1.76 1.65 35.7  -1.066  0.5413
##  April - July     2.99 2.27 38.4   1.314  0.3963
##  May - July       4.75 2.41 39.9   1.970  0.1328
## 
## tmt = No Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -2.00 2.71 35.0  -0.738  0.4654
##  April - July   nonEst   NA   NA      NA      NA
##  May - July     nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes
# stat diffs by tmt
(emmeans(tmt_diffs_mass, pairwise ~ tmt | month)$contrasts)
## month = April:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt     3.12 2.79 44.7   1.116  0.5097
##  Water Tmt - No Tmt      11.25 3.90 44.1   2.885  0.0163
##  Sham Tmt - No Tmt        8.13 3.85 43.5   2.115  0.0985
## 
## month = May:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt     6.40 3.01 52.9   2.127  0.0940
##  Water Tmt - No Tmt      14.30 4.02 47.7   3.554  0.0025
##  Sham Tmt - No Tmt        7.89 3.88 44.7   2.033  0.1160
## 
## month = July:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt     2.65 3.45 63.4   0.769  0.4448
##  Water Tmt - No Tmt     nonEst   NA   NA      NA      NA
##  Sham Tmt - No Tmt      nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes

Body Condition

tmt_diffs_SMI <- lmerTest::lmer(data = dat_repeats, SMI ~ month*tmt +
                                        (1|individual_ID))
anova(tmt_diffs_SMI, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month     282.941 141.470     2 43.652 10.7878 0.0001559 ***
## tmt        42.122  21.061     2 33.543  1.6068 0.2155784    
## month:tmt  45.205  15.068     3 41.791  1.1490 0.3405994    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
(emmeans(tmt_diffs_SMI, pairwise ~ month | tmt)$emmean)
## tmt = Water Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   34.7 1.11 64.2     32.4     36.9
##  May     39.6 1.34 65.9     37.0     42.3
##  July    31.2 1.42 66.0     28.4     34.1
## 
## tmt = Sham Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   34.2 1.04 63.4     32.2     36.3
##  May     36.2 1.11 64.8     33.9     38.4
##  July    31.5 1.64 65.3     28.2     34.8
## 
## tmt = No Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   33.4 1.79 63.4     29.8     37.0
##  May     34.4 1.79 63.4     30.8     38.0
##  July  nonEst   NA   NA       NA       NA
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
# stat diffs by month
(emmeans(tmt_diffs_SMI, pairwise ~ month | tmt)$contrasts)
## tmt = Water Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -4.99 1.62 43.0  -3.077  0.0100
##  April - July     3.42 1.69 44.8   2.022  0.1186
##  May - July       8.41 1.89 54.6   4.439  0.0001
## 
## tmt = Sham Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -1.94 1.38 37.3  -1.400  0.3514
##  April - July     2.70 1.84 46.7   1.471  0.3142
##  May - July       4.64 1.91 51.4   2.430  0.0481
## 
## tmt = No Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     -1.03 2.29 35.5  -0.449  0.6564
##  April - July   nonEst   NA   NA      NA      NA
##  May - July     nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes
# stat diffs by tmt
(emmeans(tmt_diffs_SMI, pairwise ~ tmt | month)$contrasts)
## month = April:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt    0.429 1.52 63.9   0.282  0.9572
##  Water Tmt - No Tmt      1.265 2.11 63.6   0.599  0.8211
##  Sham Tmt - No Tmt       0.837 2.07 63.4   0.404  0.9142
## 
## month = May:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt    3.479 1.74 65.6   1.997  0.1211
##  Water Tmt - No Tmt      5.224 2.24 64.6   2.332  0.0584
##  Sham Tmt - No Tmt       1.745 2.11 63.8   0.826  0.6881
## 
## month = July:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt   -0.291 2.17 65.8  -0.134  0.8937
##  Water Tmt - No Tmt     nonEst   NA   NA      NA      NA
##  Sham Tmt - No Tmt      nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes

Osmolality

tmt_diffs_osml <- lmerTest::lmer(data = dat_repeats, 
                                      osmolality_mmol_kg_mean ~ month*tmt +
                                        (1|individual_ID))
anova(tmt_diffs_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##           Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## month     4511.6 2255.81     2 43.087  4.6155 0.01527 *
## tmt        343.7  171.85     2 33.790  0.3518 0.70600  
## month:tmt  823.2  274.41     3 41.048  0.5614 0.64352  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
(emmeans(tmt_diffs_osml, pairwise ~ month | tmt)$emmean)
## tmt = Water Tmt:
##  month emmean    SE   df lower.CL upper.CL
##  April    356  7.18 62.4      342      371
##  May      379  8.29 63.8      363      396
##  July     345  8.80 64.0      327      362
## 
## tmt = Sham Tmt:
##  month emmean    SE   df lower.CL upper.CL
##  April    364  6.42 60.8      351      377
##  May      375  6.90 62.5      361      389
##  July     359 10.14 63.2      339      379
## 
## tmt = No Tmt:
##  month emmean    SE   df lower.CL upper.CL
##  April    368 12.45 62.6      343      392
##  May      376 11.13 60.8      354      399
##  July  nonEst    NA   NA       NA       NA
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
# stat diffs by month
(emmeans(tmt_diffs_osml, pairwise ~ month | tmt)$contrasts)
## tmt = Water Tmt:
##  contrast     estimate    SE   df t.ratio p.value
##  April - May    -23.01 10.04 39.7  -2.291  0.0688
##  April - July    11.58 10.62 45.6   1.090  0.5249
##  May - July      34.59 11.65 52.8   2.969  0.0123
## 
## tmt = Sham Tmt:
##  contrast     estimate    SE   df t.ratio p.value
##  April - May    -10.87  8.45 35.5  -1.286  0.4121
##  April - July     4.80 11.25 44.2   0.427  0.9048
##  May - July      15.67 11.71 48.6   1.338  0.3812
## 
## tmt = No Tmt:
##  contrast     estimate    SE   df t.ratio p.value
##  April - May     -8.74 15.05 37.3  -0.580  0.5652
##  April - July   nonEst    NA   NA      NA      NA
##  May - July     nonEst    NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes
# stat diffs by tmt
(emmeans(tmt_diffs_osml, pairwise ~ tmt | month)$contrasts)
## month = April:
##  contrast             estimate    SE   df t.ratio p.value
##  Water Tmt - Sham Tmt    -7.53  9.64 61.8  -0.781  0.7159
##  Water Tmt - No Tmt     -11.11 14.38 62.5  -0.773  0.7208
##  Sham Tmt - No Tmt       -3.58 14.01 62.3  -0.256  0.9646
## 
## month = May:
##  contrast             estimate    SE   df t.ratio p.value
##  Water Tmt - Sham Tmt     4.61 10.79 63.4   0.427  0.9045
##  Water Tmt - No Tmt       3.16 13.88 62.2   0.227  0.9719
##  Sham Tmt - No Tmt       -1.45 13.09 61.4  -0.111  0.9933
## 
## month = July:
##  contrast             estimate    SE   df t.ratio p.value
##  Water Tmt - Sham Tmt   -14.31 13.42 63.8  -1.066  0.2904
##  Water Tmt - No Tmt     nonEst    NA   NA      NA      NA
##  Sham Tmt - No Tmt      nonEst    NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes

CEWL

tmt_diffs_CEWL <- lmerTest::lmer(data = dat_repeats, 
                                      CEWL_g_m2h ~ month*tmt +
                                        (1|individual_ID))
anova(tmt_diffs_CEWL, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month     211.726 105.863     2 42.321  9.5390 0.0003805 ***
## tmt        39.478  19.739     2 34.064  1.7790 0.1841270    
## month:tmt  43.448  14.483     3 40.735  1.3049 0.2858705    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# check covariates
tmt_diffs_CEWL_check_covaries <- lmerTest::lmer(data = dat_repeats, 
                                                     CEWL_g_m2h ~ month*tmt
                                         + msmt_temp_C + msmt_VPD_kPa +
                                        (1|individual_ID))
anova(tmt_diffs_CEWL_check_covaries, type = "3", ddf = "Kenward-Roger") # not important
## Type III Analysis of Variance Table with Kenward-Roger's method
##               Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## month        158.396  79.198     2 46.720  6.8694 0.002425 **
## tmt           35.565  17.783     2 33.062  1.5441 0.228516   
## msmt_temp_C    5.959   5.959     1 61.410  0.5174 0.474666   
## msmt_VPD_kPa   7.641   7.641     1 61.763  0.6635 0.418456   
## month:tmt     42.580  14.193     3 39.151  1.2321 0.310989   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
(emmeans(tmt_diffs_CEWL, pairwise ~ month | tmt)$emmean)
## tmt = Water Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April   9.32 1.09 60.8     7.13    11.50
##  May     6.11 1.30 64.6     3.51     8.71
##  July   11.60 1.38 64.9     8.84    14.35
## 
## tmt = Sham Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April  10.04 1.02 59.2     8.00    12.08
##  May     9.30 1.09 61.9     7.12    11.48
##  July   12.98 1.58 64.1     9.83    16.13
## 
## tmt = No Tmt:
##  month emmean   SE   df lower.CL upper.CL
##  April  13.43 1.97 62.3     9.50    17.37
##  May     7.87 1.77 59.2     4.34    11.41
##  July  nonEst   NA   NA       NA       NA
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
# stat diffs by month
(emmeans(tmt_diffs_CEWL, pairwise ~ month | tmt)$contrasts)
## tmt = Water Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     3.208 1.51 40.6   2.131  0.0961
##  April - July   -2.278 1.57 42.2  -1.447  0.3264
##  May - July     -5.487 1.78 50.5  -3.075  0.0093
## 
## tmt = Sham Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     0.741 1.28 35.9   0.581  0.8312
##  April - July   -2.942 1.71 43.4  -1.719  0.2095
##  May - July     -3.683 1.79 47.3  -2.060  0.1093
## 
## tmt = No Tmt:
##  contrast     estimate   SE   df t.ratio p.value
##  April - May     5.558 2.28 37.9   2.439  0.0195
##  April - July   nonEst   NA   NA      NA      NA
##  May - July     nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes
# stat diffs by sex
(emmeans(tmt_diffs_CEWL, pairwise ~ tmt | month)$contrasts)
## month = April:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt   -0.721 1.50 60.1  -0.482  0.8800
##  Water Tmt - No Tmt     -4.112 2.25 62.0  -1.826  0.1696
##  Sham Tmt - No Tmt      -3.391 2.22 61.7  -1.529  0.2845
## 
## month = May:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt   -3.188 1.70 63.8  -1.876  0.1538
##  Water Tmt - No Tmt     -1.763 2.20 61.7  -0.803  0.7027
##  Sham Tmt - No Tmt       1.426 2.08 60.0   0.686  0.7723
## 
## month = July:
##  contrast             estimate   SE   df t.ratio p.value
##  Water Tmt - Sham Tmt   -1.385 2.10 64.8  -0.661  0.5110
##  Water Tmt - No Tmt     nonEst   NA   NA      NA      NA
##  Sham Tmt - No Tmt      nonEst   NA   NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for varying family sizes

Hydric physiology was not significantly different between treatments at any time.

Was change overall (across tmts) different from zero?

t.test(dat_tmt_change_May_only$mass_change_from_April)
## 
##  One Sample t-test
## 
## data:  dat_tmt_change_May_only$mass_change_from_April
## t = 3.828, df = 20, p-value = 0.001051
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.259063 4.274270
## sample estimates:
## mean of x 
##  2.766667
t.test(dat_tmt_change_May_only$CEWL_change_from_April)
## 
##  One Sample t-test
## 
## data:  dat_tmt_change_May_only$CEWL_change_from_April
## t = -1.8247, df = 20, p-value = 0.08303
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -4.805881  0.321135
## sample estimates:
## mean of x 
## -2.242373
t.test(dat_tmt_change_May_only$osml_change_from_April)
## 
##  One Sample t-test
## 
## data:  dat_tmt_change_May_only$osml_change_from_April
## t = 1.8095, df = 20, p-value = 0.08542
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.027158 28.566840
## sample estimates:
## mean of x 
##  13.26984

On average, lizards who were measured in both April and May increased mass, but did not experience a change in osmolality or CEWL.

Plot Delta ~ Time*Sex

Mass

ggplot(dat_tmt_change) +
  aes(x = month,
      y = mass_change_from_April,
      group = individual_ID,
      color = sex_M_F,
      shape = sex_M_F
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'Lizard Mass (g)')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes, name = "Sex") +
  savs_ggplot_theme -> delta_mass_sex
delta_mass_sex

Osmolality

ggplot(dat_tmt_change) +
  aes(x = month,
      y = osml_change_from_April,
      group = individual_ID,
      color = sex_M_F,
      shape = sex_M_F
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'Osmolality (mmol '*kg^-1*')')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes, name = "Sex") +
  savs_ggplot_theme -> delta_osml_sex
delta_osml_sex

CEWL

ggplot(dat_tmt_change) +
  aes(x = month,
      y = CEWL_change_from_April,
      group = individual_ID,
      color = sex_M_F,
      shape = sex_M_F
      ) +
  xlab("Measurement Month") + 
  ylab(expression(Delta ~ 'CEWL (g '*m^-2*' '*h^-1*')')) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1.5) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes, name = "Sex") +
  savs_ggplot_theme -> delta_CEWL_sex
delta_CEWL_sex

Export Multi-Panel Fig

Save plots:

ggarrange(delta_CEWL_sex,
          delta_osml_sex, 
          delta_mass_sex,
          ncol = 3, nrow = 1,
          #ncol = 1, nrow = 3,
          labels = c("A", "B", "C"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> delta_multi_fig_sex

ggsave(filename = "delta_multi_fig_sex.pdf",
       plot = delta_multi_fig_sex,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 200, height = 80
       #width = 80, height = 210
       )

Time*Sex Marginal Means

I want to know whether the average values for the lizards we measured changes over time. To test this, I can run a linear model, then use emmeans to get the marginal means and confidence intervals for each time point.

Do this only for lizards with repeat measurements.

info for using emmeans with interactions: https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html

Check n’s:

dat_repeats %>%
  group_by(month, sex_M_F) %>%
  summarise(n = n())
## # A tibble: 6 × 3
## # Groups:   month [3]
##   month sex_M_F     n
##   <fct> <fct>   <int>
## 1 April Female     11
## 2 April Male       22
## 3 May   Female     10
## 4 May   Male       17
## 5 July  Female      2
## 6 July  Male       12

Mass

time_sex_diffs_mass <- lmerTest::lmer(data = dat_repeats, mass_g ~ month*sex_M_F +
                                        (1|individual_ID))
anova(time_sex_diffs_mass, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month         314.003 157.002     2 38.639  9.6139 0.0004086 ***
## sex_M_F         0.032   0.032     1 38.343  0.0020 0.9647658    
## month:sex_M_F  92.617  46.309     2 38.639  2.8357 0.0709404 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
emmeans_mass <- data.frame(emmeans(time_sex_diffs_mass, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "Mass")

# stat diffs by month
time_emmeans_mass <- data.frame(emmeans(time_sex_diffs_mass, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "Mass")

# stat diffs by sex
sex_emmeans_mass <- data.frame(emmeans(time_sex_diffs_mass, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "Mass")

Body Condition

time_sex_diffs_SMI <- lmerTest::lmer(data = dat_repeats, SMI ~ month*sex_M_F +
                                        (1|individual_ID))
anova(time_sex_diffs_SMI, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month         327.19 163.597     2 45.863 15.2615 8.308e-06 ***
## sex_M_F         0.63   0.628     1 46.566  0.0587  0.809671    
## month:sex_M_F 134.56  67.280     2 45.863  6.2764  0.003896 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
emmeans_SMI <- data.frame(emmeans(time_sex_diffs_SMI, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "SMI")

# stat diffs by month
time_emmeans_SMI <- data.frame(emmeans(time_sex_diffs_SMI, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "SMI")

# stat diffs by sex
sex_emmeans_SMI <- data.frame(emmeans(time_sex_diffs_SMI, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "SMI")

Hematocrit

time_sex_diffs_hct <- lmerTest::lmer(data = dat_repeats, 
                                      hematocrit_percent ~ month*sex_M_F +
                                        (1|individual_ID))
anova(time_sex_diffs_hct, type = "2", ddf = "Kenward-Roger")
## Type II Analysis of Variance Table with Kenward-Roger's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month         817.60  408.80     2 44.421 15.8405 6.412e-06 ***
## sex_M_F       169.86  169.86     1 31.941  6.5855   0.01518 *  
## month:sex_M_F  21.74   10.87     2 45.442  0.4211   0.65885    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
emmeans_hct <- data.frame(emmeans(time_sex_diffs_hct, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "Hematocrit")

# stat diffs by month
time_emmeans_hct <- data.frame(emmeans(time_sex_diffs_hct, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "Hematocrit")

# stat diffs by sex
sex_emmeans_hct <- data.frame(emmeans(time_sex_diffs_hct, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "Hematocrit")

Osmolality

time_sex_diffs_osml <- lmerTest::lmer(data = dat_repeats, 
                                      osmolality_mmol_kg_mean ~ month*sex_M_F +
                                        (1|individual_ID))
anova(time_sex_diffs_osml, type = "2", ddf = "Kenward-Roger")
## Type II Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## month          8073.9  4037.0     2 48.038  9.6102 0.0003086 ***
## sex_M_F       11642.2 11642.2     1 31.800 27.7401 9.302e-06 ***
## month:sex_M_F  1301.2   650.6     2 49.692  1.5476 0.2228392    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
emmeans_osml <- data.frame(emmeans(time_sex_diffs_osml, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "Osmolality")

# stat diffs by month
time_emmeans_osml <- data.frame(emmeans(time_sex_diffs_osml, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "Osmolality")

# stat diffs by sex
sex_emmeans_osml <- data.frame(emmeans(time_sex_diffs_osml, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "Osmolality")

CEWL

time_sex_diffs_CEWL <- lmerTest::lmer(data = dat_repeats, 
                                      CEWL_g_m2h ~ month*sex_M_F +
                                        (1|individual_ID))
anova(time_sex_diffs_CEWL, type = "2", ddf = "Kenward-Roger")
## Type II Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## month         180.393  90.197     2 44.819  8.0478 0.001032 **
## sex_M_F       103.055 103.055     1 32.419  9.1999 0.004737 **
## month:sex_M_F  18.921   9.461     2 46.365  0.8438 0.436564   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# check covariates
time_sex_diffs_CEWL_check_covaries <- lmerTest::lmer(data = dat_repeats, 
                                                     CEWL_g_m2h ~ month*sex_M_F
                                         + msmt_temp_C + msmt_VPD_kPa +
                                        (1|individual_ID))
anova(time_sex_diffs_CEWL_check_covaries, type = "2", ddf = "Kenward-Roger") # not important at all
## Type II Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)   
## month         101.714  50.857     2 50.265  4.3866 0.01754 * 
## sex_M_F        89.529  89.529     1 34.970  7.7316 0.00868 **
## msmt_temp_C     0.133   0.133     1 64.720  0.0115 0.91496   
## msmt_VPD_kPa    0.090   0.090     1 64.849  0.0078 0.93004   
## month:sex_M_F  17.518   8.759     2 46.384  0.7555 0.47547   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
59/0.09 # 655.5
## [1] 655.5556
54/0.14 # 385.7
## [1] 385.7143
# mean values
emmeans_CEWL <- data.frame(emmeans(time_sex_diffs_CEWL, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "CEWL")

# stat diffs by month
time_emmeans_CEWL <- data.frame(emmeans(time_sex_diffs_CEWL, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "CEWL")

# stat diffs by sex
sex_emmeans_CEWL <- data.frame(emmeans(time_sex_diffs_CEWL, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "CEWL")

Check Assumptions

plot(time_sex_diffs_CEWL)

plot(time_sex_diffs_osml)

plot(time_sex_diffs_mass)

plot(time_sex_diffs_SMI)

Merge dfs

emmeans_all <- emmeans_CEWL %>%
  rbind(emmeans_osml) %>%
  rbind(emmeans_mass) %>%
  rbind(emmeans_SMI) %>%
  mutate(month_sex = factor(paste(month, sex_M_F),
                            levels = c("April Female", "April Male",
                                       "May Female", "May Male",
                                       "July Female", "July Male")))

# all pairwise stats together
emmean_time_stats <- time_emmeans_CEWL %>%
  rbind(time_emmeans_osml) %>%
  rbind(time_emmeans_mass) %>%
  rbind(time_emmeans_SMI)
emmean_sex_stats <- sex_emmeans_CEWL %>%
  rbind(sex_emmeans_osml) %>%
  rbind(sex_emmeans_mass) %>%
  rbind(sex_emmeans_SMI)

Plot MMeans ~ Time*Sex

These are plots of the marginal means, only for the dataset of lizards with repeat measurements.

Mass

ggplot() +
  geom_jitter(data = dat_repeats,
             aes(x = month_sex,
                 y = mass_g,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = emmeans_all[emmeans_all$measurement == "Mass",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = emmeans_all[emmeans_all$measurement == "Mass",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Lizard Mass (g)") +
  # sex differences = NONE
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 61, label = "A", size = 3) +
  annotate(geom = "text", x = 4, y = 61, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 61, label = "A", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 61, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 61, label = "b", size = 3) +
  annotate(geom = "text", x = 5, y = 61, label = "c", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels) +
  ylim(17, 61) +
  savs_ggplot_theme -> mass_emmean_plot
mass_emmean_plot

Body Condition

ggplot() +
  geom_jitter(data = dat_repeats,
             aes(x = month_sex,
                 y = SMI,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = emmeans_all[emmeans_all$measurement == "SMI",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = emmeans_all[emmeans_all$measurement == "SMI",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Body Condition (g')") +
  # sex differences
  annotate(geom = "text", x = 3.5, y = 45, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 54, label = "AB", size = 3) +
  annotate(geom = "text", x = 4, y = 54, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 54, label = "B", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 54, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 54, label = "b", size = 3) +
  annotate(geom = "text", x = 5, y = 54, label = "a", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels) +
  ylim(20, 55) +
  savs_ggplot_theme -> SMI_emmean_plot
SMI_emmean_plot

Osmolality

ggplot() +
  geom_jitter(data = dat_repeats,
             aes(x = month_sex,
                 y = osmolality_mmol_kg_mean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = emmeans_all[emmeans_all$measurement == "Osmolality",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = emmeans_all[emmeans_all$measurement == "Osmolality",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
  # sex differences
  annotate(geom = "text", x = 1.5, y = 416, label = "*", size = 6) +
  annotate(geom = "text", x = 3.5, y = 416, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 430, label = "A", size = 3) +
  annotate(geom = "text", x = 4, y = 430, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 430, label = "B", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 430, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 430, label = "a", size = 3) +
  annotate(geom = "text", x = 5, y = 430, label = "a", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels) +
  scale_y_continuous(limits = c(320, 440), breaks = c(320, 360, 400, 440)) +
  savs_ggplot_theme -> osml_emmean_plot
osml_emmean_plot

CEWL

ggplot() +
  geom_jitter(data = dat_repeats,
             aes(x = month_sex,
                 y = CEWL_g_m2h,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
             position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = emmeans_all[emmeans_all$measurement == "CEWL",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = emmeans_all[emmeans_all$measurement == "CEWL",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                    color = sex_M_F),
                width = .4) +
  # sex differences
  annotate(geom = "text", x = 3.5, y = 16, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 19, label = "A", size = 3) +
  annotate(geom = "text", x = 4, y = 19, label = "B", size = 3) +
  annotate(geom = "text", x = 6, y = 19, label = "A", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 19, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 19, label = "a", size = 3) +
  annotate(geom = "text", x = 5, y = 19, label = "a", size = 3) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels) +
  savs_ggplot_theme -> CEWL_emmean_plot
CEWL_emmean_plot

This is only for lizards that were recaptured and measured more than once, so I think these changes may reflect true seasonality in their water balance.

Export Repeats-Only Multi-Panel Fig

Put in one plot together:

ggarrange(CEWL_emmean_plot, 
          osml_emmean_plot, 
          mass_emmean_plot,
          SMI_emmean_plot,
          ncol = 4, nrow = 1,
#          ncol = 1, nrow = 3,
          labels = c("A", "B", "C", "D"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> szn_change_multi_fig
#szn_change_multi_fig

ggsave(filename = "szn_change_multi_fig.pdf",
       plot = szn_change_multi_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 200, height = 80
       #width = 80, height = 210
       )

FULL Time*Sex Marginal Means

I’m going to do literally the same exact thing as the stats and figures in the section above, but now for all of my data, not just for the data on lizards with repeat measurements.

Check n’s:

dat %>%
  group_by(month, sex_M_F) %>%
  summarise(n = n())
## # A tibble: 6 × 3
## # Groups:   month [3]
##   month sex_M_F     n
##   <fct> <fct>   <int>
## 1 April Female     26
## 2 April Male       41
## 3 May   Female     15
## 4 May   Male       22
## 5 July  Female      3
## 6 July  Male       12

Mass

full_time_sex_diffs_mass <- lmerTest::lmer(data = dat, mass_g ~ month*sex_M_F +
                                        (1|individual_ID))
full_time_sex_diffs_mass_aov <- data.frame(anova(full_time_sex_diffs_mass, 
                                                 type = "2", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))

# mean values
full_emmeans_mass <- data.frame(emmeans(full_time_sex_diffs_mass, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "Mass")

# stat diffs by month
full_time_emmeans_mass <- data.frame(emmeans(full_time_sex_diffs_mass, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "Mass")

# stat diffs by sex
full_sex_emmeans_mass <- data.frame(emmeans(full_time_sex_diffs_mass, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "Mass")

Body Condition

full_time_sex_diffs_SMI <- lmerTest::lmer(data = dat, SMI ~ month*sex_M_F +
                                        (1|individual_ID))
full_time_sex_diffs_SMI_aov <- data.frame(anova(full_time_sex_diffs_SMI, 
                                          type = "3", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))

# double check that type 2 vs 3 anovas are essentially the same results
full_time_sex_diffs_SMI_aov2 <- data.frame(anova(full_time_sex_diffs_SMI, 
                                          type = "2", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))
full_time_sex_diffs_SMI_aov
##                  Sum.Sq   Mean.Sq NumDF     DenDF   F.value       Pr..F.
## month         400.48274 200.24137     2  79.56813 18.583496 2.384786e-07
## sex_M_F        11.62646  11.62646     1 102.97466  1.079132 3.013248e-01
## month:sex_M_F 138.38877  69.19439     2  79.56813  6.421618 2.597665e-03
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
full_time_sex_diffs_SMI_aov2
##                   Sum.Sq    Mean.Sq NumDF    DenDF   F.value       Pr..F.
## month         348.170763 174.085381     2 71.15328 16.157819 1.639477e-06
## sex_M_F         8.594231   8.594231     1 74.72959  0.797690 3.746521e-01
## month:sex_M_F 138.388775  69.194387     2 79.56813  6.421618 2.597665e-03
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
# mean values
full_emmeans_SMI <- data.frame(emmeans(full_time_sex_diffs_SMI, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "SMI")

# stat diffs by month
full_time_emmeans_SMI <- data.frame(emmeans(full_time_sex_diffs_SMI, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "SMI")

# stat diffs by sex
full_sex_emmeans_SMI <- data.frame(emmeans(full_time_sex_diffs_SMI, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "SMI")

Osmolality

full_time_sex_diffs_osml <- lmerTest::lmer(data = dat, 
                                      osmolality_mmol_kg_mean ~ month*sex_M_F +
                                        (1|individual_ID))
full_time_sex_diffs_osml_aov <- data.frame(anova(full_time_sex_diffs_osml, 
                                           type = "3", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))
full_time_sex_diffs_osml_aov2 <- data.frame(anova(full_time_sex_diffs_osml, 
                                           type = "2", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))
full_time_sex_diffs_osml_aov
##                  Sum.Sq   Mean.Sq NumDF    DenDF   F.value      Pr..F.
## month         6748.7057 3374.3528     2 91.86749 5.8321054 0.004125822
## sex_M_F       1441.1831 1441.1831     1 88.14015 2.4917205 0.118028949
## month:sex_M_F  770.2048  385.1024     2 91.86749 0.6655966 0.516427511
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
full_time_sex_diffs_osml_aov2
##                  Sum.Sq   Mean.Sq NumDF    DenDF   F.value       Pr..F.
## month         9627.0309 4813.5154     2 87.00381 8.3213727 0.0004933249
## sex_M_F       5545.4928 5545.4928     1 69.65080 9.5878298 0.0028217467
## month:sex_M_F  770.2048  385.1024     2 91.86749 0.6655966 0.5164275115
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
# mean values
full_emmeans_osml <- data.frame(emmeans(full_time_sex_diffs_osml, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "Osmolality")

# stat diffs by month
full_time_emmeans_osml <- data.frame(emmeans(full_time_sex_diffs_osml, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "Osmolality")

# stat diffs by sex
full_sex_emmeans_osml <- data.frame(emmeans(full_time_sex_diffs_osml, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "Osmolality")

CEWL

full_time_sex_diffs_CEWL <- lmerTest::lmer(data = dat, 
                                      CEWL_g_m2h ~ month*sex_M_F +
                                        (1|individual_ID))
full_time_sex_diffs_CEWL_aov <- data.frame(anova(full_time_sex_diffs_CEWL, 
                                           type = "3", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))
full_time_sex_diffs_CEWL_aov2 <- data.frame(anova(full_time_sex_diffs_CEWL, 
                                           type = "2", ddf = "Kenward-Roger")) %>%
  mutate(explanatory = factor(row.names(.)))
full_time_sex_diffs_CEWL_aov
##                  Sum.Sq  Mean.Sq NumDF    DenDF  F.value      Pr..F.
## month         119.80851 59.90425     2 85.07588 5.625373 0.005075837
## sex_M_F        34.27453 34.27453     1 98.76587 3.219121 0.075842334
## month:sex_M_F  46.88928 23.44464     2 85.07588 2.201594 0.116890977
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
full_time_sex_diffs_CEWL_aov2
##                  Sum.Sq  Mean.Sq NumDF    DenDF  F.value       Pr..F.
## month         181.78030 90.89015     2 77.41997 8.536304 0.0004465354
## sex_M_F        55.19526 55.19526     1 73.45087 5.184032 0.0257090750
## month:sex_M_F  46.88928 23.44464     2 85.07588 2.201594 0.1168909770
##                 explanatory
## month                 month
## sex_M_F             sex_M_F
## month:sex_M_F month:sex_M_F
# check covariates
full_time_sex_diffs_CEWL_check_covaries <- lmerTest::lmer(data = dat, 
                                                     CEWL_g_m2h ~ month*sex_M_F
                                         + msmt_temp_C + msmt_VPD_kPa +
                                        (1|individual_ID))
anova(full_time_sex_diffs_CEWL_check_covaries, 
      type = "2", ddf = "Kenward-Roger")
## Type II Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## month         185.755  92.878     2  89.095  8.5769 0.0003921 ***
## sex_M_F        46.450  46.450     1  74.944  4.2908 0.0417603 *  
## msmt_temp_C     6.098   6.098     1 109.982  0.5633 0.4545481    
## msmt_VPD_kPa    9.009   9.009     1 109.977  0.8322 0.3636216    
## month:sex_M_F  32.730  16.365     2  86.160  1.5113 0.2264115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean values
full_emmeans_CEWL <- data.frame(emmeans(full_time_sex_diffs_CEWL, 
                                   pairwise ~ month | sex_M_F)$emmean) %>%
  mutate(measurement = "CEWL")

# stat diffs by month
full_time_emmeans_CEWL <- data.frame(emmeans(full_time_sex_diffs_CEWL, 
                                   pairwise ~ month | sex_M_F)$contrasts) %>%
  mutate(measurement = "CEWL")

# stat diffs by sex
full_sex_emmeans_CEWL <- data.frame(emmeans(full_time_sex_diffs_CEWL, 
                                   pairwise ~ sex_M_F | month)$contrasts) %>%
  mutate(measurement = "CEWL")

Check Assumptions

plot(full_time_sex_diffs_CEWL)

plot(full_time_sex_diffs_osml)

plot(full_time_sex_diffs_mass)

plot(full_time_sex_diffs_SMI)

Merge dfs

full_emmeans_all <- full_emmeans_CEWL %>%
  rbind(full_emmeans_osml) %>%
  rbind(full_emmeans_mass) %>%
  rbind(full_emmeans_SMI) %>%
  mutate(month_sex = factor(paste(month, sex_M_F),
                            levels = c("April Female", "April Male",
                                       "May Female", "May Male",
                                       "July Female", "July Male")))
full_emmeans_F <- full_emmeans_all %>%
  dplyr::filter(sex_M_F == "Female")
full_emmeans_M <- full_emmeans_all %>%
  dplyr::filter(sex_M_F == "Male")

# all pairwise stats together
full_emmean_time_stats <- full_time_emmeans_CEWL %>%
  rbind(full_time_emmeans_osml) %>%
  rbind(full_time_emmeans_mass) %>%
  rbind(full_time_emmeans_SMI)
full_emmean_sex_stats <- full_sex_emmeans_CEWL %>%
  rbind(full_sex_emmeans_osml) %>%
  rbind(full_sex_emmeans_mass) %>%
  rbind(full_sex_emmeans_SMI)

# anova stats - type 2 is better for looking at interactions
full_time_sex_diffs_all <- (full_time_sex_diffs_CEWL_aov2) %>%
  rbind((full_time_sex_diffs_osml_aov2)) %>%
  #rbind((full_time_sex_diffs_mass_aov2)) %>%
  rbind((full_time_sex_diffs_SMI_aov2)) %>%
  mutate(response = factor(c(rep("CEWL", 3), rep("Plasma Osmolality", 3),
                      #rep("Mass", 3), 
                      rep("Body Condition", 3))))
summary(full_time_sex_diffs_all)
##      Sum.Sq            Mean.Sq             NumDF           DenDF      
##  Min.   :   8.594   Min.   :   8.594   Min.   :1.000   Min.   :69.65  
##  1st Qu.:  55.195   1st Qu.:  55.195   1st Qu.:1.000   1st Qu.:73.45  
##  Median : 181.780   Median :  90.890   Median :2.000   Median :77.42  
##  Mean   :1857.972   Mean   :1240.613   Mean   :1.667   Mean   :78.88  
##  3rd Qu.: 770.205   3rd Qu.: 385.102   3rd Qu.:2.000   3rd Qu.:85.08  
##  Max.   :9627.031   Max.   :5545.493   Max.   :2.000   Max.   :91.87  
##     F.value            Pr..F.                 explanatory              response
##  Min.   : 0.6656   Min.   :0.0000016   month        :3    Body Condition   :3  
##  1st Qu.: 2.2016   1st Qu.:0.0004933   month:sex_M_F:3    CEWL             :3  
##  Median : 6.4216   Median :0.0028217   sex_M_F      :3    Plasma Osmolality:3  
##  Mean   : 6.4304   Mean   :0.1155601                                           
##  3rd Qu.: 8.5363   3rd Qu.:0.1168910                                           
##  Max.   :16.1578   Max.   :0.5164275

Export Model Results

full_time_sex_diffs_4exp <- full_time_sex_diffs_all %>%
  mutate(partial_sum_of_squares = round(Sum.Sq, 0), 
         F_statistic = round(F.value, 2),
         DenDF = round(DenDF, 0),
         F_full = paste(F_statistic, " (", NumDF, ",", DenDF, ")", sep = ""),
         p_value = round(Pr..F., 4)) %>%
  dplyr::select(response, explanatory,
                partial_sum_of_squares,
                #NumDF, DenDF,
                #F_statistic, 
                F_full,
                p_value)
write.csv(full_time_sex_diffs_4exp, 
          "./results_statistics/hydric_Fstats_by_monthsex.csv")

FULL Plot MMeans Time*Sex

Mass

ggplot() +
  geom_jitter(data = dat,
             aes(x = month_sex,
                 y = mass_g,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_all[full_emmeans_all$measurement == "Mass",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_all[full_emmeans_all$measurement == "Mass",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Lizard Mass (g)") +
  # sex differences
  annotate(geom = "text", x = 5.5, y = 54, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 61, label = "AB", size = 3) +
  annotate(geom = "text", x = 4, y = 61, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 61, label = "B", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 61, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 61, label = "b", size = 3) +
  annotate(geom = "text", x = 5, y = 61, label = "a", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels_full) +
  ylim(19, 61) +
  savs_ggplot_theme -> full_mass_emmean_plot
full_mass_emmean_plot

Mass - F

ggplot() +
  geom_jitter(data = dat_F,
             aes(x = month,
                 y = mass_g,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_F[full_emmeans_F$measurement == "Mass",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_F[full_emmeans_F$measurement == "Mass",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Lizard Mass (g)") +
  # sex differences
  annotate(geom = "text", x = 3, y = 54, label = "*", size = 6) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 61, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 61, label = "b", size = 3) +
  annotate(geom = "text", x = 3, y = 61, label = "a", size = 3) +
  scale_color_manual(values = Fem, name = "Sex") +
  scale_fill_manual(values = Fem, name = "Sex") +
  scale_shape_manual(values = 19, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_F) +
  ylim(19, 61) +
  savs_ggplot_theme -> full_mass_emmean_plot_F
full_mass_emmean_plot_F

ggsave(filename = "szn_change_mass_F_standalone.pdf",
       plot = full_mass_emmean_plot_F,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Mass - M

ggplot() +
  geom_jitter(data = dat_M,
             aes(x = month,
                 y = mass_g,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_M[full_emmeans_M$measurement == "Mass",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_M[full_emmeans_M$measurement == "Mass",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Lizard Mass (g)") +
  # sex differences
  annotate(geom = "text", x = 3, y = 54, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 1, y = 61, label = "a,b", size = 3) +
  annotate(geom = "text", x = 2, y = 61, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 61, label = "b", size = 3) +
  scale_color_manual(values = Male, name = "Sex") +
  scale_fill_manual(values = Male, name = "Sex") +
  scale_shape_manual(values = 15, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_M) +
  ylim(19, 61) +
  savs_ggplot_theme -> full_mass_emmean_plot_M
full_mass_emmean_plot_M

ggsave(filename = "szn_change_mass_M_standalone.pdf",
       plot = full_mass_emmean_plot_M,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Body Condition

ggplot() +
  geom_jitter(data = dat,
             aes(x = month_sex,
                 y = SMI,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_all[full_emmeans_all$measurement == "SMI",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_all[full_emmeans_all$measurement == "SMI",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Body Condition (g')") +
  # sex differences
  annotate(geom = "text", x = 1.5, y = 44, label = "*", size = 6) +
  annotate(geom = "text", x = 3.5, y = 44, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 49, label = "A", size = 3) +
  annotate(geom = "text", x = 4, y = 49, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 49, label = "B", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 49, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 49, label = "b", size = 3) +
  annotate(geom = "text", x = 5, y = 49, label = "a", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels_full) +
  ylim(19, 50) +
  savs_ggplot_theme -> full_SMI_emmean_plot
full_SMI_emmean_plot

Body Cond - F

ggplot() +
  geom_jitter(data = dat_F,
             aes(x = month,
                 y = SMI,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_F[full_emmeans_F$measurement == "SMI",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_F[full_emmeans_F$measurement == "SMI",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Body Condition (g')") +
  # sex differences
  annotate(geom = "text", x = 1, y = 47, label = "*", size = 6) +
  annotate(geom = "text", x = 2, y = 47, label = "*", size = 6) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 50, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 50, label = "b", size = 3) +
  annotate(geom = "text", x = 3, y = 50, label = "a", size = 3) +
  scale_color_manual(values = Fem, name = NULL) +
  scale_fill_manual(values = Fem, name = NULL) +
  scale_shape_manual(values = 21, name = NULL) +
  scale_x_discrete(labels = my_labels_full_F) +
  ylim(19, 50) +
  savs_ggplot_theme +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 10, unit = "pt"),
        legend.margin = margin(20,0,0,0, unit = "pt")) -> full_SMI_emmean_plot_F
full_SMI_emmean_plot_F

ggsave(filename = "szn_change_bodycond_F_standalone.pdf",
       plot = full_SMI_emmean_plot_F,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Body Cond - M

ggplot() +
  geom_jitter(data = dat_M,
             aes(x = month,
                 y = SMI,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_M[full_emmeans_M$measurement == "SMI",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_M[full_emmeans_M$measurement == "SMI",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab("Body Condition (g')") +
  # sex differences
  annotate(geom = "text", x = 1, y = 47, label = "*", size = 6) +
  annotate(geom = "text", x = 2, y = 47, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 1, y = 50, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 50, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 50, label = "b", size = 3) +
  scale_color_manual(values = Male, name = NULL) +
  scale_fill_manual(values = Male, name = NULL) +
  scale_shape_manual(values = 22, name = NULL) +
  scale_x_discrete(labels = my_labels_full_M) +
  ylim(19, 50) +
  savs_ggplot_theme  +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"),
        legend.margin = margin(20,0,0,0, unit = "pt")) -> full_SMI_emmean_plot_M
full_SMI_emmean_plot_M

ggsave(filename = "szn_change_bodycond_M_standalone.pdf",
       plot = full_SMI_emmean_plot_M,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Osmolality

ggplot() +
  geom_jitter(data = dat,
             aes(x = month_sex,
                 y = osmolality_mmol_kg_mean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_all[full_emmeans_all$measurement == "Osmolality",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_all[full_emmeans_all$measurement == "Osmolality",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
  # sex differences
  annotate(geom = "text", x = 1.5, y = 422, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 445, label = "AB", size = 3) +
  annotate(geom = "text", x = 4, y = 445, label = "A", size = 3) +
  annotate(geom = "text", x = 6, y = 445, label = "B", size = 3) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 445, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 445, label = "a", size = 3) +
  annotate(geom = "text", x = 5, y = 445, label = "a", size = 3) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels_full) +
  scale_y_continuous(limits = c(320, 450), breaks = c(320, 360, 400, 440)) +
  savs_ggplot_theme -> full_osml_emmean_plot
full_osml_emmean_plot

Osmolality - F

ggplot() +
  geom_jitter(data = dat_F,
             aes(x = month,
                 y = osmolality_mmol_kg_mean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_F[full_emmeans_F$measurement == "Osmolality",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_F[full_emmeans_F$measurement == "Osmolality",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
  # sex differences
  annotate(geom = "text", x = 1, y = 435, label = "*", size = 6) +
  # time differences for FEMALES
  annotate(geom = "text", x = 1, y = 448, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 448, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 448, label = "a", size = 3) +
  scale_color_manual(values = Fem, name = "Sex") +
  scale_fill_manual(values = Fem, name = "Sex") +
  scale_shape_manual(values = 21, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_F) +
  scale_y_continuous(limits = c(320, 450), breaks = c(320, 360, 400, 440)) +
  savs_ggplot_theme +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"),
        legend.margin = margin(0,0,0,0, unit = "pt")) -> full_osml_emmean_plot_F
full_osml_emmean_plot_F

ggsave(filename = "szn_change_osml_F_standalone.pdf",
       plot = full_osml_emmean_plot_F,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Osmolality - M

ggplot() +
  geom_jitter(data = dat_M,
             aes(x = month,
                 y = osmolality_mmol_kg_mean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_M[full_emmeans_M$measurement == "Osmolality",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_M[full_emmeans_M$measurement == "Osmolality",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                 color = sex_M_F),
                width = .4) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
  # sex differences
  annotate(geom = "text", x = 1, y = 435, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 1, y = 448, label = "a,b", size = 3) +
  annotate(geom = "text", x = 2, y = 448, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 448, label = "b", size = 3) +
  scale_color_manual(values = Male, name = "Sex") +
  scale_fill_manual(values = Male, name = "Sex") +
  scale_shape_manual(values = 22, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_M) +
  scale_y_continuous(limits = c(320, 450), breaks = c(320, 360, 400, 440)) +
  savs_ggplot_theme +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"),
        legend.margin = margin(0,0,0,0, unit = "pt")) -> full_osml_emmean_plot_M
full_osml_emmean_plot_M

ggsave(filename = "szn_change_osml_M_standalone.pdf",
       plot = full_osml_emmean_plot_M,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

CEWL

ggplot() +
  geom_jitter(data = dat,
             aes(x = month_sex,
                 y = CEWL_g_m2h,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
             position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_all[full_emmeans_all$measurement == "CEWL",],
             aes(x = month_sex,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_all[full_emmeans_all$measurement == "CEWL",],
                aes(x = month_sex,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                    color = sex_M_F),
                width = .4) +
  # sex differences
  annotate(geom = "text", x = 3.5, y = 20, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 2, y = 24, label = "A", size = 3) +
  annotate(geom = "text", x = 4, y = 24, label = "B", size = 3) +
  annotate(geom = "text", x = 6, y = 24, label = "A", size = 3) +
  # time differences for FEMALES = none
  annotate(geom = "text", x = 1, y = 24, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 24, label = "a", size = 3) +
  annotate(geom = "text", x = 5, y = 24, label = "a", size = 3) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
  ylim(0, 25) +
  scale_color_manual(values = my_colors, name = "Sex") +
  scale_fill_manual(values = my_colors, name = "Sex") +
  scale_shape_manual(values = my_shapes_open, name = "Sex") +
  scale_x_discrete(labels = my_labels_full) +
  savs_ggplot_theme -> full_CEWL_emmean_plot
full_CEWL_emmean_plot

CEWL - F

ggplot() +
  geom_jitter(data = dat_F,
             aes(x = month,
                 y = CEWL_g_m2h,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
             position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_F[full_emmeans_F$measurement == "CEWL",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_F[full_emmeans_F$measurement == "CEWL",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                    color = sex_M_F),
                width = .4) +
  # sex differences
  annotate(geom = "text", x = 2, y = 22.5, label = "*", size = 6) +
  # time differences for FEMALES = none
  annotate(geom = "text", x = 1, y = 25, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 25, label = "a", size = 3) +
  annotate(geom = "text", x = 3, y = 25, label = "a", size = 3) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
  ylim(0, 25) +
  scale_color_manual(values = Fem, name = "Sex") +
  scale_fill_manual(values = Fem, name = "Sex") +
  scale_shape_manual(values = 21, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_F) +
  savs_ggplot_theme +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 4.5, unit = "pt"),
        legend.margin = margin(0,0,0,0, unit = "pt")) -> full_CEWL_emmean_plot_F
full_CEWL_emmean_plot_F

ggsave(filename = "szn_change_CEWL_F_standalone.pdf",
       plot = full_CEWL_emmean_plot_F,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

CEWL - M

ggplot() +
  geom_jitter(data = dat_M,
             aes(x = month,
                 y = CEWL_g_m2h,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 0.3,
             size = 0.9,
             position = position_jitter(height = 0, width = 0.2)) +
  geom_point(data = full_emmeans_M[full_emmeans_M$measurement == "CEWL",],
             aes(x = month,
                 y = emmean,
                 shape = sex_M_F,
                 color = sex_M_F,
                 fill = sex_M_F),
             alpha = 1,
             size = 2) +
  geom_errorbar(data = full_emmeans_M[full_emmeans_M$measurement == "CEWL",],
                aes(x = month,
                    y = emmean, 
                    ymin = lower.CL,
                    ymax = upper.CL,
                    color = sex_M_F),
                width = .4) +
  # sex differences
  annotate(geom = "text", x = 2, y = 22.5, label = "*", size = 6) +
  # time differences for MALES
  annotate(geom = "text", x = 1, y = 25, label = "a", size = 3) +
  annotate(geom = "text", x = 2, y = 25, label = "b", size = 3) +
  annotate(geom = "text", x = 3, y = 25, label = "a", size = 3) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
  ylim(0, 25) +
  scale_color_manual(values = Male, name = "Sex") +
  scale_fill_manual(values = Male, name = "Sex") +
  scale_shape_manual(values = 22, name = "Sex") +
  scale_x_discrete(labels = my_labels_full_M) +
  savs_ggplot_theme +
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"),
        legend.margin = margin(0,0,0,0, unit = "pt")) -> full_CEWL_emmean_plot_M
full_CEWL_emmean_plot_M

ggsave(filename = "szn_change_CEWL_M_standalone.pdf",
       plot = full_CEWL_emmean_plot_M,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 60, height = 80
       )

Arrange and Export M&F Separate Panel Plot

# CEWL
full_CEWL_emmean_plot_M_multi <- full_CEWL_emmean_plot_M + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  labs(y = NULL)
full_CEWL_emmean_plot_F_multi <- full_CEWL_emmean_plot_F + 
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none")
# osmolality
full_osml_emmean_plot_M_multi <- full_osml_emmean_plot_M + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  labs(y = NULL)
full_osml_emmean_plot_F_multi <- full_osml_emmean_plot_F + 
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none")
# body condition
full_SMI_emmean_plot_M_multi <- full_SMI_emmean_plot_M + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  labs(y = NULL)

# arrange M & F together
ggarrange(full_CEWL_emmean_plot_F_multi, full_CEWL_emmean_plot_M_multi,
          full_osml_emmean_plot_F_multi, full_osml_emmean_plot_M_multi,
          full_SMI_emmean_plot_F, full_SMI_emmean_plot_M_multi,
          ncol = 2, nrow = 3,
          widths = c(1.36, 1), # make x-axis widths match when exported
          heights = c(1, 1, 1.36), # heights
          labels = c("(a)", "", "(b)", "", "(c)", ""),
          label.y = 1.02,
          label.x = -0.05,
          font.label = list(color = "black", face = "bold", size = 10)
          ) -> full_szn_change_multi_fig_best

ggsave(filename = "szn_change_multi_fig_official.pdf",
       plot = full_szn_change_multi_fig_best,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 220 # don't change bc might mess up relative sizing
       )

Export Full Data Multi-Panel Fig

ggarrange(full_CEWL_emmean_plot, 
          full_osml_emmean_plot, 
          #full_mass_emmean_plot,
          full_SMI_emmean_plot,
          ncol = 3, nrow = 1,
          #ncol = 1, nrow = 3,
          labels = c("A", "B", "C"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> full_szn_change_multi_fig_h
#full_szn_change_multi_fig

ggsave(filename = "szn_change_multi_fig_full_horizontal.pdf",
       plot = full_szn_change_multi_fig_h,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 220, height = 80
       #width = 80, height = 210
       )

Export Full & Repeats Multi-Panel Fig

ggarrange(full_CEWL_emmean_plot, 
          full_osml_emmean_plot, 
          full_mass_emmean_plot,
          full_SMI_emmean_plot,
          CEWL_emmean_plot, 
          osml_emmean_plot, 
          mass_emmean_plot,
          SMI_emmean_plot,
          ncol = 4, nrow = 2,
#          ncol = 1, nrow = 3,
          labels = c("full", "full", "full", "full",
                     "rep", "rep", "rep", "rep"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> full_vs_rep_szn_change_multi_fig
#full_vs_rep_szn_change_multi_fig

#ggsave(filename = "szn_change_multi_fig_full_vs_rep.pdf",
 #      plot = full_vs_rep_szn_change_multi_fig,
  #     path = "./results_figures",
   #    device = "pdf",
    #   dpi = 600,
     #  units = "mm",
      # width = 210, height = 160
       ##width = 80, height = 210
       #)

Reproductive Effort

Is reproductive effort related to water balance?

Egg Investment

summary(dat_repeat_females)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 11:39:00.00   Min.   :2021-04-23   Length:18          
##  1st Qu.:2021-04-23 15:07:30.00   1st Qu.:2021-04-23   Class :character   
##  Median :2021-04-24 12:51:30.00   Median :2021-05-01   Mode  :character   
##  Mean   :2021-04-30 06:59:56.25   Mean   :2021-04-30                      
##  3rd Qu.:2021-05-08 09:32:30.00   3rd Qu.:2021-05-08                      
##  Max.   :2021-05-08 12:27:00.00   Max.   :2021-05-08                      
##  NA's   :2                                                                
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:18          F-01   :2     Female:18   Min.   :31.00   Min.   : 90.0  
##  Class :character   F-03   :2     Male  : 0   1st Qu.:35.60   1st Qu.: 95.0  
##  Mode  :character   F-05   :2                 Median :40.00   Median :100.0  
##                     F-06   :2                 Mean   :41.73   Mean   :100.7  
##                     F-11   :2                 3rd Qu.:46.75   3rd Qu.:106.5  
##                     F-12   :2                 Max.   :58.80   Max.   :114.0  
##                     (Other):6                                                
##         tmt     tmt_mass_change_g capture_to_msmt   hematocrit_percent
##  Water Tmt: 8   Min.   :-3.0000   Min.   :  25.00   Min.   :23.00     
##  Sham Tmt :10   1st Qu.:-1.4000   1st Qu.:  55.25   1st Qu.:30.25     
##  No Tmt   : 0   Median :-0.1000   Median : 147.00   Median :32.50     
##                 Mean   :-0.5941   Mean   : 225.88   Mean   :33.78     
##                 3rd Qu.: 0.0000   3rd Qu.: 247.00   3rd Qu.:37.00     
##                 Max.   : 1.8000   Max.   :1638.00   Max.   :44.00     
##                 NA's   :1         NA's   :2                           
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :28.00   Min.   : 1.493   Min.   :23.13   Min.   :12.20  
##  1st Qu.:33.25   1st Qu.: 8.986   1st Qu.:30.47   1st Qu.:14.19  
##  Median :34.00   Median :10.320   Median :31.59   Median :14.90  
##  Mean   :33.72   Mean   :11.109   Mean   :30.79   Mean   :16.29  
##  3rd Qu.:35.00   3rd Qu.:13.559   3rd Qu.:31.89   3rd Qu.:18.02  
##  Max.   :36.00   Max.   :21.673   Max.   :32.23   Max.   :28.17  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month   week  
##  Min.   :2.830   Min.   :2.033   Min.   :329.5           April:9   1 :9  
##  1st Qu.:4.358   1st Qu.:3.692   1st Qu.:338.0           May  :9   3 :9  
##  Median :4.646   Median :3.916   Median :345.0           July :0   12:0  
##  Mean   :4.463   Mean   :3.751   Mean   :350.9                           
##  3rd Qu.:4.725   3rd Qu.:4.054   3rd Qu.:362.9                           
##  Max.   :4.818   Max.   :4.167   Max.   :390.0                           
##                                                                          
##     week_num ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1   Min.   :2021-04-24   Not Gravid: 6   Min.   :3.000  
##  1st Qu.:1   1st Qu.:2021-04-24   Gravid    :12   1st Qu.:3.750  
##  Median :2   Median :2021-05-01                   Median :4.000  
##  Mean   :2   Mean   :2021-05-01                   Mean   :3.833  
##  3rd Qu.:3   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3   Max.   :2021-05-08                   Max.   :5.000  
##                                                   NA's   :6      
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage  
##  Min.   :0.640           Min.   :1.820                small round:1  
##  1st Qu.:0.865           1st Qu.:2.542                large round:6  
##  Median :1.095           Median :3.185                soft       :2  
##  Mean   :1.173           Mean   :3.292                soft oblong:1  
##  3rd Qu.:1.488           3rd Qu.:4.040                firm oblong:1  
##  Max.   :1.800           Max.   :4.740                hard oblong:1  
##  NA's   :6               NA's   :6                    NA's       :6  
##    dev_point         SMI               month_sex date_chunk_pretty
##  Min.   :1.00   Min.   :29.24   April Female:9   Apr 23-25:9      
##  1st Qu.:2.00   1st Qu.:34.63   April Male  :0   May 7-9  :9      
##  Median :2.00   Median :36.73   May Female  :9   Jul 14   :0      
##  Mean   :2.75   Mean   :37.83   May Male    :0                    
##  3rd Qu.:3.25   3rd Qu.:41.85   July Female :0                    
##  Max.   :5.00   Max.   :45.63   July Male   :0                    
##  NA's   :6                                                        
##  percent_water_mass_change n     
##  Min.   :-6.834            2:18  
##  1st Qu.:-4.242            1: 0  
##  Median :-0.250                  
##  Mean   :-1.539                  
##  3rd Qu.: 0.000                  
##  Max.   : 3.462                  
##  NA's   :1
dat_repeat_females_eggs <- dat_repeat_females %>% 
  dplyr::filter(gravid_Y_N == "Gravid") %>% 
  dplyr::filter(month == "May") # so we don't have any repeats

summary(lm(data = dat_repeat_females_eggs, number_eggs ~ mass_g))
## 
## Call:
## lm(formula = number_eggs ~ mass_g, data = dat_repeat_females_eggs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2855 -0.2647 -0.1891  0.1467  0.7541 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.11895    0.69935   1.600   0.1607   
## mass_g       0.05975    0.01554   3.846   0.0085 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4103 on 6 degrees of freedom
## Multiple R-squared:  0.7114, Adjusted R-squared:  0.6633 
## F-statistic: 14.79 on 1 and 6 DF,  p-value: 0.008501
summary(lm(data = dat_repeat_females_eggs, number_eggs ~ SVL_mm))
## 
## Call:
## lm(formula = number_eggs ~ SVL_mm, data = dat_repeat_females_eggs)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38064 -0.30282 -0.00921  0.28214  0.44492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -5.43434    2.07306  -2.621  0.03951 * 
## SVL_mm       0.09173    0.02066   4.439  0.00438 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.369 on 6 degrees of freedom
## Multiple R-squared:  0.7666, Adjusted R-squared:  0.7277 
## F-statistic: 19.71 on 1 and 6 DF,  p-value: 0.00438
summary(lm(data = dat_repeat_females_eggs, number_eggs ~ SMI))
## 
## Call:
## lm(formula = number_eggs ~ SMI, data = dat_repeat_females_eggs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9073 -0.3127 -0.0824  0.2596  0.9036 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.04386    2.40577   0.018    0.986
## SMI          0.09238    0.05970   1.548    0.173
## 
## Residual standard error: 0.6457 on 6 degrees of freedom
## Multiple R-squared:  0.2853, Adjusted R-squared:  0.1661 
## F-statistic: 2.395 on 1 and 6 DF,  p-value: 0.1727
anova(lm(data = dat_repeat_females_eggs, number_eggs ~ mass_g))
## Analysis of Variance Table
## 
## Response: number_eggs
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## mass_g     1 2.4899 2.48991   14.79 0.008501 **
## Residuals  6 1.0101 0.16835                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(data = dat_repeat_females_eggs, number_eggs ~ SVL_mm))
## Analysis of Variance Table
## 
## Response: number_eggs
##           Df  Sum Sq Mean Sq F value  Pr(>F)   
## SVL_mm     1 2.68307 2.68307  19.706 0.00438 **
## Residuals  6 0.81693 0.13616                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Repeat Individuals

I have 9 females measured at both April and May time points. Let’s see if their physiology changes consistently based on becoming gravid.

First, calculate the change in CEWL, osmolality, and mass for each of those females between April to May.

fem_starting_values <- dat_repeat_females %>%
  # get the values for April measurements only 
  dplyr::filter(month == "April") %>%
  dplyr::select(individual_ID,
                CEWL_April = CEWL_g_m2h,
                osml_April = osmolality_mmol_kg_mean,
                mass_April = mass_g,
                SMI_April = SMI)
fem_change <- dat_repeat_females %>%
  # add starting values
  left_join(fem_starting_values, by = 'individual_ID') %>%
  # make sure we have complete data
  dplyr::filter(complete.cases(mass_April, CEWL_April)) %>%
  # calculate difference for each msmt versus in April
  mutate(mass_change_from_April = mass_g - mass_April,
         SMI_change_from_April = SMI - SMI_April,
         CEWL_change_from_April = CEWL_g_m2h - CEWL_April,
         osml_change_from_April = osmolality_mmol_kg_mean - osml_April)
summary(fem_change)
##  capture_date_time                 capture_date        radio_collar_serial
##  Min.   :2021-04-23 11:39:00.00   Min.   :2021-04-23   Length:18          
##  1st Qu.:2021-04-23 15:07:30.00   1st Qu.:2021-04-23   Class :character   
##  Median :2021-04-24 12:51:30.00   Median :2021-05-01   Mode  :character   
##  Mean   :2021-04-30 06:59:56.25   Mean   :2021-04-30                      
##  3rd Qu.:2021-05-08 09:32:30.00   3rd Qu.:2021-05-08                      
##  Max.   :2021-05-08 12:27:00.00   Max.   :2021-05-08                      
##  NA's   :2                                                                
##   PIT_tag_ID        individual_ID   sex_M_F       mass_g          SVL_mm     
##  Length:18          F-01   :2     Female:18   Min.   :31.00   Min.   : 90.0  
##  Class :character   F-03   :2     Male  : 0   1st Qu.:35.60   1st Qu.: 95.0  
##  Mode  :character   F-05   :2                 Median :40.00   Median :100.0  
##                     F-06   :2                 Mean   :41.73   Mean   :100.7  
##                     F-11   :2                 3rd Qu.:46.75   3rd Qu.:106.5  
##                     F-12   :2                 Max.   :58.80   Max.   :114.0  
##                     (Other):6                                                
##         tmt     tmt_mass_change_g capture_to_msmt   hematocrit_percent
##  Water Tmt: 8   Min.   :-3.0000   Min.   :  25.00   Min.   :23.00     
##  Sham Tmt :10   1st Qu.:-1.4000   1st Qu.:  55.25   1st Qu.:30.25     
##  No Tmt   : 0   Median :-0.1000   Median : 147.00   Median :32.50     
##                 Mean   :-0.5941   Mean   : 225.88   Mean   :33.78     
##                 3rd Qu.: 0.0000   3rd Qu.: 247.00   3rd Qu.:37.00     
##                 Max.   : 1.8000   Max.   :1638.00   Max.   :44.00     
##                 NA's   :1         NA's   :2                           
##  cloacal_temp_C    CEWL_g_m2h      msmt_temp_C    msmt_RH_percent
##  Min.   :28.00   Min.   : 1.493   Min.   :23.13   Min.   :12.20  
##  1st Qu.:33.25   1st Qu.: 8.986   1st Qu.:30.47   1st Qu.:14.19  
##  Median :34.00   Median :10.320   Median :31.59   Median :14.90  
##  Mean   :33.72   Mean   :11.109   Mean   :30.79   Mean   :16.29  
##  3rd Qu.:35.00   3rd Qu.:13.559   3rd Qu.:31.89   3rd Qu.:18.02  
##  Max.   :36.00   Max.   :21.673   Max.   :32.23   Max.   :28.17  
##                                                                  
##     e_s_kPa       msmt_VPD_kPa   osmolality_mmol_kg_mean   month   week  
##  Min.   :2.830   Min.   :2.033   Min.   :329.5           April:9   1 :9  
##  1st Qu.:4.358   1st Qu.:3.692   1st Qu.:338.0           May  :9   3 :9  
##  Median :4.646   Median :3.916   Median :345.0           July :0   12:0  
##  Mean   :4.463   Mean   :3.751   Mean   :350.9                           
##  3rd Qu.:4.725   3rd Qu.:4.054   3rd Qu.:362.9                           
##  Max.   :4.818   Max.   :4.167   Max.   :390.0                           
##                                                                          
##     week_num ultrasound_date           gravid_Y_N  number_eggs   
##  Min.   :1   Min.   :2021-04-24   Not Gravid: 6   Min.   :3.000  
##  1st Qu.:1   1st Qu.:2021-04-24   Gravid    :12   1st Qu.:3.750  
##  Median :2   Median :2021-05-01                   Median :4.000  
##  Mean   :2   Mean   :2021-05-01                   Mean   :3.833  
##  3rd Qu.:3   3rd Qu.:2021-05-08                   3rd Qu.:4.000  
##  Max.   :3   Max.   :2021-05-08                   Max.   :5.000  
##                                                   NA's   :6      
##  largest_egg_diameter_cm largest_egg_circumference_cm         stage  
##  Min.   :0.640           Min.   :1.820                small round:1  
##  1st Qu.:0.865           1st Qu.:2.542                large round:6  
##  Median :1.095           Median :3.185                soft       :2  
##  Mean   :1.173           Mean   :3.292                soft oblong:1  
##  3rd Qu.:1.488           3rd Qu.:4.040                firm oblong:1  
##  Max.   :1.800           Max.   :4.740                hard oblong:1  
##  NA's   :6               NA's   :6                    NA's       :6  
##    dev_point         SMI               month_sex date_chunk_pretty
##  Min.   :1.00   Min.   :29.24   April Female:9   Apr 23-25:9      
##  1st Qu.:2.00   1st Qu.:34.63   April Male  :0   May 7-9  :9      
##  Median :2.00   Median :36.73   May Female  :9   Jul 14   :0      
##  Mean   :2.75   Mean   :37.83   May Male    :0                    
##  3rd Qu.:3.25   3rd Qu.:41.85   July Female :0                    
##  Max.   :5.00   Max.   :45.63   July Male   :0                    
##  NA's   :6                                                        
##  percent_water_mass_change n        CEWL_April       osml_April   
##  Min.   :-6.834            2:18   Min.   : 1.493   Min.   :329.5  
##  1st Qu.:-4.242            1: 0   1st Qu.: 8.940   1st Qu.:333.3  
##  Median :-0.250                   Median :10.353   Median :338.0  
##  Mean   :-1.539                   Mean   :11.304   Mean   :344.8  
##  3rd Qu.: 0.000                   3rd Qu.:13.640   3rd Qu.:347.7  
##  Max.   : 3.462                   Max.   :21.673   Max.   :390.0  
##  NA's   :1                                                        
##    mass_April      SMI_April     mass_change_from_April SMI_change_from_April
##  Min.   :31.00   Min.   :29.24   Min.   :-2.100         Min.   :-0.6421      
##  1st Qu.:33.00   1st Qu.:32.19   1st Qu.: 0.000         1st Qu.: 0.0000      
##  Median :40.00   Median :35.13   Median : 0.000         Median : 0.0000      
##  Mean   :39.44   Mean   :34.93   Mean   : 2.289         Mean   : 2.8985      
##  3rd Qu.:44.00   3rd Qu.:37.02   3rd Qu.: 4.425         3rd Qu.: 3.8356      
##  Max.   :52.00   Max.   :41.30   Max.   : 9.000         Max.   :14.2653      
##                                                                              
##  CEWL_change_from_April osml_change_from_April
##  Min.   :-8.2493        Min.   :-50.667       
##  1st Qu.:-1.5265        1st Qu.:  0.000       
##  Median : 0.0000        Median :  0.000       
##  Mean   :-0.1947        Mean   :  6.111       
##  3rd Qu.: 0.0000        3rd Qu.: 11.625       
##  Max.   : 7.8247        Max.   : 45.000       
## 

Plot Indiv Change

ggplot(fem_change) +
  aes(x = month,
      y = mass_change_from_April,
      group = individual_ID,
      color = gravid_Y_N
      ) +
  geom_point() +
  geom_line() +
  theme_bw()

ggplot(fem_change) +
  aes(x = month,
      y = CEWL_change_from_April,
      group = individual_ID,
      color = gravid_Y_N
      ) +
  geom_point() +
  geom_line() +
  theme_bw()

ggplot(fem_change) +
  aes(x = month,
      y = osml_change_from_April,
      group = individual_ID,
      color = gravid_Y_N
      ) +
  geom_point() +
  geom_line() +
  theme_bw()

boxplot(fem_change$osml_change_from_April)

boxplot(fem_change$mass_change_from_April)

On average, mass and osmolality both increase. There’s only one lizard who was not gravid at both time points, so we can’t compare change based on becoming gravid or not. I think it is reasonable to conclude that on average, gravid females gained mass and became less hydrated. But, we cannot say that this was DUE to gravidity, because we have nothing to compare it to

Stats Indiv Change

Run stats tests on the change rates for only the females who were gravid in May, and use only the actual delta values, in their rows for May data.

# check osml outlier range
fem_change %>%
  group_by(month) %>%
  summarise(mean(osml_change_from_April), sd(osml_change_from_April))
## # A tibble: 2 × 3
##   month `mean(osml_change_from_April)` `sd(osml_change_from_April)`
##   <fct>                          <dbl>                        <dbl>
## 1 April                            0                            0  
## 2 May                             12.2                         29.0
# for stats
fem_change_stats <- fem_change %>%
  # don't use the one female still non-gravid in May
  dplyr::filter(individual_ID != "F-12") %>%
  # only use the >resulting< change values
  dplyr::filter(month == "May") %>%
  # remove the crazy low outlier based on sd calculated above
  dplyr::filter(osml_change_from_April > (12-(29*2)))

ttests:

t.test(fem_change_stats$CEWL_change_from_April)
## 
##  One Sample t-test
## 
## data:  fem_change_stats$CEWL_change_from_April
## t = -0.59229, df = 6, p-value = 0.5753
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -5.873479  3.584194
## sample estimates:
## mean of x 
## -1.144643
t.test(fem_change_stats$osml_change_from_April)
## 
##  One Sample t-test
## 
## data:  fem_change_stats$osml_change_from_April
## t = 2.8997, df = 6, p-value = 0.02734
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   3.28688 38.80836
## sample estimates:
## mean of x 
##  21.04762
t.test(fem_change_stats$mass_change_from_April)
## 
##  One Sample t-test
## 
## data:  fem_change_stats$mass_change_from_April
## t = 3.5281, df = 6, p-value = 0.0124
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.514744 8.370970
## sample estimates:
## mean of x 
##  4.942857
t.test(fem_change_stats$SMI_change_from_April)
## 
##  One Sample t-test
## 
## data:  fem_change_stats$SMI_change_from_April
## t = 2.7067, df = 6, p-value = 0.03526
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   0.5695776 11.2980848
## sample estimates:
## mean of x 
##  5.933831

Stats Likelihood

Does anything correlate with the likelihood of being/becoming gravid?

Start with trying each possible logistic model to predict gravidity:

repro_logit1 <- lme4::glmer(data = dat_repro_females, 
                            gravid_Y_N ~ mass_g +
                            (1|individual_ID),
                            family = "binomial")

repro_logit2 <- lme4::glmer(data = dat_repro_females, 
                            gravid_Y_N ~ SMI +
                            (1|individual_ID),
                            family = "binomial")

repro_logit3 <- lme4::glmer(data = dat_repro_females, 
                            gravid_Y_N ~ osmolality_mmol_kg_mean +
                            (1|individual_ID),
                            family = "binomial")

repro_logit4 <- lme4::glmer(data = dat_repro_females, 
                            gravid_Y_N ~ CEWL_g_m2h +
                            (1|individual_ID),
                            family = "binomial")


car::Anova(repro_logit1, type = 2, 
           test.statistic = "Chisq", error.estimate = "pearson")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: gravid_Y_N
##         Chisq Df Pr(>Chisq)
## mass_g 2.3158  1     0.1281
car::Anova(repro_logit2, type = 2, 
           test.statistic = "Chisq", error.estimate = "pearson")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: gravid_Y_N
##      Chisq Df Pr(>Chisq)
## SMI 1.3948  1     0.2376
car::Anova(repro_logit3, type = 2, 
           test.statistic = "Chisq", error.estimate = "pearson")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: gravid_Y_N
##                          Chisq Df Pr(>Chisq)
## osmolality_mmol_kg_mean 0.2219  1     0.6376
car::Anova(repro_logit4, type = 2, 
           test.statistic = "Chisq", error.estimate = "pearson")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: gravid_Y_N
##             Chisq Df Pr(>Chisq)
## CEWL_g_m2h 0.0045  1     0.9468

What about an ANOVA?

Plot Likelihood Gravid ~ Mass

dat_repro_females %>%
  # dummy-fy gravid variable
  mutate(gravid_Y_N = factor(gravid_Y_N, levels = c("Not Gravid", "Gravid"),
                             labels = c(0,1)),
         # goes 1-2, so fix to be 0-1
         gravid_Y_N = as.numeric(gravid_Y_N)-1) %>%
  # nooooowwwww I can plot and the stat_smooth will actually work
  ggplot() +
  aes(x = mass_g, 
      y = gravid_Y_N) +
  geom_jitter(alpha=.5, 
              position = position_jitter(height = 0.1, width = 0)
              ) +
  # must make response variable 0-1 in order for this to work
  stat_smooth(geom = "smooth",
              method = "glm", 
              formula = y~x,
              se = FALSE, 
              method.args = list(family = "quasibinomial")
              ) +
  ylab("") +
  scale_y_continuous(labels = c("Not Gravid", "", "", "", "Gravid"), 
                     limits = c(0,1)) +
  xlab("Body Mass (g)") +
  savs_ggplot_theme

Summary Figure

dodger <- position_dodge(width = 0.2) # make dodge amt for both points and lines

dat_repro_females %>%
  # only use data for April-May
  dplyr::filter(month %in% c("April", "May")) %>% 
  ggplot() +
    aes(x = month,
        y = mass_g,
        color = gravid_Y_N,
        shape = gravid_Y_N,
        alpha = gravid_Y_N,
        group = individual_ID) +
  geom_line(position = dodger) +
  geom_point(size = 2,
              position = dodger, #position_jitter(height = 0.01, width = 0),
             # alpha = 0.5
             ) +
  labs(y = "Body Mass (g)",
       x = NULL) +
  scale_color_manual(values = c(Fem_nongrav, Fem_grav), name = NULL) +
  scale_shape_manual(values = c(19, 19), name = NULL) +
  scale_alpha_manual(values = c(0.3, 0.7), name = NULL) +
  savs_ggplot_theme +
  theme(legend.position = "right",
        plot.margin = margin(t = 1.5, r = 0.5, b = 0.5, l = 0.5, unit = "pt"),
        legend.margin = margin(0,0,0,0, unit = "pt")) -> when_gravid
when_gravid

ggsave(filename = "female_gravidity.pdf",
       plot = when_gravid,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       #width = 200, height = 80
       width = 80, height = 60
       )